AD-A226  652 


m  a 


FH  F.  COPY 


Observation  and  Inversion  of  Seismo-Acoustic  Waves  in  a 
Complex  Arctic  Ice  Environment 


Bruce  Edward  Miller 

D.S.,  United  States  Naval  Academy  (1975) 

Submitted  in  partial  fulfillment  of  the 
requirements  for  the  dual  degrees  of 

OCEAN  ENGINEER 


TYJir 

AS  A  k  V* 

01  HLECTrl  F 

SEP  2 13  1990$ 

-  rd 


ft/PS 


MASTER  OF  SCIENCE  IN  OCEAN  ENGINEERING 

at  the 

MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 

and  the  fj0a,  2^3*  tf-G-etto 

WOODS  HOLE  OCEANOGRAPHIC  INSTITUTION 


U 


September  1990  DfS'.T  07 /.*.  ; : : : ;  T 

©  Bruce  Edward  Miller,  1990  AppiovoJ 

i  Pi^TASv.aaa  urJuaiiud 

The  author  hereby  grants  to  M.I.T.,  W.H.O.I.  and  the  U.S.  Government 
permission  to  reproduce  and  to  distribute  copies  of  this  thesis  document  in  whole 

or  in  part.  , 


Signature  of  Author. 


Joint  Program  in  Oceanographic  Engineering 
Massachusetts  Institute  of  Technology 
Woods  Hole  Oceanographic  Institution 
~~  August  10,  1990 


Certified  by. 


certified  by. . ... .(. .  .S.,y. .-.A- 

(  )  . 


Dr.  Henrik  Schmidt 
Massachusetts  Institute  of  Technology 
Thesis  Supervisor 


Dr.  James  F.  Lynch 


— "■  ^00^8  ^olc  QjWfl,nographic  Institution 

Accepted  by  ... .  yMS-  ])r  Vy  'Kcnjall  Melville 

Chairman, Joint  Committee  for  Oceanographic  Engineering 
Massachusetts  Institute  of  Tcchnology/Woods  Hole  Oceanographic  Institution 


Observation  and  Inversion  of  Seismo-Acoustic  Waves  in  a  Complex 

Arctic  Ice  Environment 
by 

Bruce  Edward  Miller 

Submitted  to  the  Massachusetts  Institute  of  Technology/ 

Woods  Hole  Oceanographic  Institution 
Joint  Program  in  Oceanographic  Engineering 
on  August  10,  1990,  in  partial  fulfillment  of  the 
requirements  for  the  degrees  of 
Ocean  Engineer 
and 

Master  of  Science  in  Ocean  Engineering 


Abstract 

The  propagation  of  low  frequency  seismo-acoustic  waves  in  the  Arctic  Ocean  ice 
canopy  is  examined  through  the  analysis  of  hydrophone  and  geophone  data  sets  collected 
in  1987  at  an  ice  camp  designated  PRUDEX  in  the  Beaufort  Sea. 

Study  of  the  geophone  time  series  generated  by  under-ice  explosive  detonations  reveals 
not  only  the  expected  longitudinal  and  flexural  waves  in  the  ice  plate,  but  also  an  unex¬ 
pected  horizontally-polarized  transverse  (SH)  wave  arriving  at  a  higher  amplitude  than 
the  other  wave  types.  The  travel  paths  of  all  three  observed  wave  types  are  found  to  be 
refracted  in  the  horizontal  plane  along  a  line  coincident  with  a  known  ridge  separating 
the  ice  canopy  locally  into  two  distinct  half-plates,  the  first  of  thin  first  year  ice  and  the 
second  of  thicker  multi-year  ice.  The  origin  of  the  SII  wave  appears  to  be  near  the  detona¬ 
tion  and  not  associated  with  the  interaction  of  longitudinal,  flexural  or  waterborne  waves 
with  the  ridge  line.  The  need  to  determine  the  exact  location  of  each  detonation  from  the 
received  time  series  highlights  the  dramatic  superiority  of  geophones  over  hydrophones  in 
this  application,  as  does  the  ability  to  detect  the  anomalous  SH  waves  and  the  refracted 
ray  paths,  neither  of  which^are  visible  in  the  hydrophone  data. 

Inversion  of  the  geophone  data  sets  for  the  low  frequency  elastic  parameters  of  the  ice 
is  conducted  initially  by  treating  the  ice  as  a  single  homogeneous  isotropic  plate  to  demon¬ 
strate  the  power  of  SAFARI  numerical  modeling  in  this  application.  A  modified  stationary 
phase  approach  is  then  used  to  extend  SAFARI  modeling  to  invert  the  data  sets  for  the 
elastic  parameters  of  the  two  ice  half-plates  simultaneously.  'The  compressional/shear 
bulk  wave  speeds  estimated  in  the  half-plates,  3500/1750  m/s  in  the  multi-year  ice  and 
3000/1590  m/s  in  the  new  ice,  are  comparable  to  previously  obtained  values;  however,  the 
compressional/shear  attenuation  values  in  the  two  half-plates,  1.0/2.99  dB/A  and  1.0/2.67 
dB/A,  respectively,  are  somewhat  greater  than  previously  measured  values  and  four  times 
greater  than  estimates  extrapolated  from  high  frequency  data. 

Thesis  Supervisor:  Dr,  Henrik  Schmidt 
Massachusetts  Institute  of  Technology 
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Chapter  1 
Introduction 


The  introduction  describes  the  motivation  and  objectives  for  this  thesis,  and 
reviews  its  organization  and  content  by  chapter. 

1.1  Motivation 

The  Arctic  Ocean  has  grown  increasingly  important  to  national  strategic  interests 
in  recent  years;  yet  our  understanding  of  the  Arctic  Ocean  environment  has  lagged  far 
behind  that  of  the  other  major  ocean  systems.  In  particular,  the  modeling  of  low 
frequency  acoustic  propagation  under  the  sea  ice  canopy  in  the  Arctic  Ocean  has  proven 
to  be  an  elusive  problem  [1].  The  difficulty  has  not  been  in  general  that  the  necessary 
tools  to  do  this  modeling  are  unavailable.  For  instance,  Schmidt’s  Fast  Field  algorithm, 
SAFARI  [2],  has  proven  to  be  a  very  capable  package  for  solving  propagation 
problems  in  a  complex  seismo-acoustic  environment  such  as  is  presented  by  the  deep 
Arctic.  The  difficulty  has  been  that  very  little  work  has  been  done  to  obtain 
measurements  of  the  starting  parameters  crucial  to  computing  this  propagation 
accurately,  i.e.,  the  elastic  parameters  of  the  ice  canopy  -  compressional  and  shear  wave 
bulk  velocities  and  attenuation  factors.  As  a  result,  previous  modeling  has  been  based 
largely  on  parameters  measured  in  the  laboratory  or  extrapolated  from  somewhat  similar 
environments  (freshwater  lake  ice,  glacial  ice,  etc.).  Recent  investigations  have 
suggested  that  parameters  so  obtained  do  not  accurately  reflect  the  Arctic  environment. 
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and  that  measurements  of  the  elastic  parameters  of  the  arctic  ice  are  needed  in  situ 
before  we  can  hope  to  model  acoustic  propagation  in  the  Arctic  Ocean  accurately. 

1.2  Thesis  Objectives 

The  overall  objective  of  this  work  is  to  improve  our  ability  to  accurately  model 
acoustic  propagation  in  the  Arctic  Ocean.  To  this  end,  data  sets  obtained  in  1987  during 
seismo-acoustic  propagation  experiments  conducted  at  an  ice  camp  in  the  Beaufort  Sea 
designated  PRUDEX  are  studied  extensively.  The  initial  intent  of  this  study  was  simply 
to  apply  advanced  modeling  techniques  to  the  problem  of  inverting  the  hydrophone  and 
geonhone  data  as  received  at  the  PRUDEX  arrays  to  obtain  accurate  measurements  of 
the  elastic  parameters.  Although  inversion  to  obtain  elastic  parameters  remains  the 
focus  of  this  work,  unexpected  phenomena  observed  in  the  propagation  data  have  served 
to  partially  frustrate  the  immedia.c  goal  of  obtaining  the  elastic  parameters  by  increasing 
the  difficulty  of  the  inversion,  while  simultaneously  contributing  to  the  overall 
understanding  of  seismo-acoustic  propagation  in  the  Arctic  by  disclosing  mechanisms 
in  the  propagation  previously  unobserved  or  unsuspected.  In  particular,  arrivals 
characteristic  of  the  refraction  of  all  types  of  propagating  waves  at  a  linear  discontinuity 
in  the  horizontal  plane  of  the  ice  plate  are  presented.  The  author  also  attempts  with 
apparent  success  to  model  this  refraction  and  extend  the  inversion  to  obtain  the  elastic 
parameters  of  two  separate  kinds  of  ice  cover,  annual  ice  and  multi-year  ice, 
simultaneously.  More  importantly,  completely  unexpected  horizontally  polarized 
transverse  (SH  mode)  waves  are  presented  in  the  propagation  data.  Although  the 
existing  theories  of  seismo-acoustic  propagation  in  an  elastic  plate  have  no  mechanism 
by  which  an  underwater  explosion  can  generate  SH  waves  in  a  sheet  of  ice  floating  over 
that  explosion,  these  waves  are  present  in  the  PRUDEX  data  sets  at  amplitudes  greater 


than  any  other  wave  type  observeu.  Incse  waves  are  important  because  their  creation 
by  out-of-plane  scattering  of  other  wave  types  may  be  a  significant  and  heretofore 
unknown  mechanism  for  the  attenuation  of  acoustic  energy  entering  the  ice  cover  from 
the  water.  In  addition  to  identifying  their  existence  in  the  PRUDEX  data,  as  a  first  step 
toward  understanding  their  origin,  the  author  investigates  possible  source  locations  for 
the  SH  waves  relative  the  explosive  source,  the  other  wave  types  in  the  ice,  and  known 
discontinuities  in  the  ice  cover. 

1.3  Thesis  Content 

Chapter  2  of  this  thesis  lays  the  foundation  for  later  work  by  reviewing  the 
theory  of  wave  propagation  in  a  thin  elastic  plate  under  various  conditions,  focusing  on 
the  development  of  the  three  wave  types  commonly  observed  in  such  plates,  the 
longitudinal  plate  wave,  the  flexural  wave,  and  the  transversely  polarized  SH  wave.  The 
second  chapter  then  reviews  briefly  the  numerical  modeling  tool  used  throughout  this 
work,  Schmidt  s  SAFARI  aigorithm  [3],  Chapter  2  concludes  with  a  short  discussion 
of  the  approach  to  the  attenuation  of  elastic  waves  employed  in  this  work. 

The  third  chapter  introduces  the  reader  to  the  experimental  data  as  obtained  at 
the  PRUDEX  ice  camp.  The  background  for  the  experiments  is  reviewed  to  make  clear 
the  need  to  determine  through  the  analysis  of  seismo-acoustic  propagation  data  some 
parameters  which  could  have  been  measured  precisely  during  the  experiment.  The 
nature  of  the  acoustic  source  pulse  used  during  the  experiments,  a  key  factor  in  later 
analysis,  is  described  and  explained  using  data  received  at  a  hydrophone  array.  Finally, 
the  occurrence  of  the  three  principal  wave  types  is  identified  in  the  data  received  at  the 
experiment’s  geopiione  array. 

Chapter  4  exists  principally  because  the  location  and  time  of  me  underwater 
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explosive  detonations  used  to  excite  seismo-acoustic  propagation  in  the  ice  cover  were 
not  recorded  relative  to  the  receiving  arrays  during  the  PRUDEX  experiments.  In  the 
process  of  determining  shot  time  and  location  for  use  in  later  inversion,  this  chapter 
brings  out  the  somewhat  startling  result  that  a  simple  array  of  four  3-axis  geophones  can 
be  a  much  more  effective  tool  for  locating  underwater  sources  than  a  larger  hydrophone 
array  in  the  water  below  the  geophones.  Geophone  data  is  used  in  Chapter  4  not  only 
to  determine  a  much  more  accurate  source  location  than  is  available  using  hydrophone 
data,  but  also  to  identify  and  analyze  the  refraction  at  the  joint  between  two  abutting  ice 
half-plates  of  all  wave  types  propagating  in  the  horizontal  plane  of  the  ice  sheet  . 
Chapter  4  also  reviews  the  evidence  available  to  help  identify  the  origin  of  the  high 
amplitude  SH  waves  which  are  visible  only  in  the  geophone  data. 

Chapter  5  begins  by  explaining  the  fundamentals  of  the  process  of  inverting 
response  data  for  the  elastic  parameters  of  the  propagating  media,  and  then  reviews 
previous  work  done  to  determine  those  parameters  in  arctic  sea  ice.  The  description  and 
results  of  the  inversion  obtained  by  treating  the  ice  canopy  as  a  single  homogeneous 
isotropic  plate  follow.  These  results  serve  to  demonstrate  the  potential  of  SAFARI 
modeling  of  wave  propagation;  although  based  on  the  work  of  Chapter  4,  the  ice  is  more 
accurately  modeled  as  two  abutting  half-plates  with  significantly  different  elastic 
parameters,  and  the  results  obtained  by  the  single  plate  model  are  of  questionable 
accuracy.  To  solve  this  problem.  Chapter  5  introduces  a  method  for  using  stationary 
phase  analysis  |4)  in  a  somewhat  modified  form  to  extend  two-dimensional  SAFARI 
to  model  propagation  of  the  flexural  wave  in  the  range  dependent  environment  of  the 
two  abutting  half-plates.  The  results  of  inversion  using  this  modeling  technique  are  then 
presented. 

Chapter  6  summarizes  the  results  of  Chapters  3,  4  and  5  and  establishes  their 
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importance  relative  to  seismo-acoustic  propagation  in  the  Arctic  Ocean.  Chapter  6  also 
comments  on  additional  work  needed  in  this  area. 
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Chapter  2 
Theory 

Prior  to  beginning  the  evaluation  and  analysis  of  any  experimental  data,  it  is 
important  to  have  a  firm  underpinning  in  the  theory  behind  that  data.  This  chapter  will 
review  briefly  the  fundamental  theory  of  propagation  of  elastic  waves  in  an  ice  plate  and 
the  equations  which  characterize  that  propagation.  It  will  then  describe  the  principal 
tool  used  in  this  paper  to  solve  these  equations  numerically,  Schmidt’s  SAFARI 
algorithms  [3]. 

2.1  Propagation  of  Elastic  Waves  in  a  Plate 

An  understanding  of  the  unique  characteristics  of  wave  propagation  in  a  thin  elastic 
plate  is  essential  to  the  study  of  seismo-acoustic  waves  in  the  arctic  ice.  Three 
fundamental  wave  types,  longitudinal  waves,  flexural  waves  and  horizontally-polarized 
transverse  (SH)  waves  are  commonly  observed  propagating  in  floating  ice  sheets 
[5]  [6].  It  is  useful  to  look  first  at  the  origins  of  the  three  wave  types  in  a  free  elastic 
plate  and  then  extend  those  results  to  a  plate  bounded  by  a  liquid  half-space  on  one  side 
and  a  vacuum  (or  air)  on  the  other. 

2.1.1  Elastic  Waves  in  a  Free  Plate 

Consider  waves  propagating  in  the  positive  x  direction  in  a  laterally  infinite 
homogeneous  isotropic  plate  bounded  on  both  sides  by  vacuum  as  shown  in  Figure  2-1. 
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Figure  2-1:  Two  dimensional  geometry  for  infinite  plate  of  thickness  2h 
bounded  by  vacuum  above  and  vacuum  or  liquid  below. 

If  displacements  in  the  plate  satisfy  the  Navier  equation, 

p u  =  /+(A+2p)V(V-m)  -  pVx(Vxw)  , 

where  X  and  p  are  the  Lame  constants,  then  Lame’s  theorem  states  [4]  that  a  scalar 
potential  4>  and  a  vector  potential  \| /  exist  which  satisfy 


u  -  V<J)  +  Vxvjf 

(2.2) 

o 

II 

> 

(2.3) 

vty  =  -U 

(2.4) 

£t2 

(2.5) 

p2 
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where 


i  A.  +2ix 

a  =  - — 

(2.6) 

P 

=  . 

(2.7) 

P 

As  a  trial  solution,  consider  plane  wave  solutions  of  the  form 

4>  = 

(2.8) 

If  =  . 

(2.9) 

Substitution  of  these  solutions  into  (2.2)  through  (2.5)  yields  the  following  set  of  four 

equations  in  eight  unknowns  for  the  potentials, 

<j)  = 

(2.10) 

(2.11) 

ljf  = 

(2.12) 

(2.13) 

The  wave  nature  of  the  potentials,  and  thus  the  motion,  is  clearly  seen  in  these 
equations. 

(2.10)  through  (2.13)  can  be  used  to  develop  seismo-acoustic  propagation  in  any 
layered  (two  dimensional)  environment,  and  various  types  of  body  and  surface  waves 
arise  depending  upon  the  media  and  the  boundary  conditions.  The  quantities  a  and  (3 
are  the  compressional  (P-wave)  and  shear  (SV/SH-waves)  bulk  velocities;  kx  is  the 
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horizontal  wavenumber  and  and  k./P  are  the  vertical  wavenumbers  for  compressional 
and  shear  waves.  No  attempt  is  made  in  this  paper  to  review  the  theory  supporting  any 
propagation  other  than  that  occurring  in  a  thin  elastic  plate;  although  references  to  other 
wave  types  are  used  to  relate  propagation  in  the  thin  elastic  plate  to  other  results  for 
readers  familiar  with  seismo-acoustics. 

From  this  point,  one  of  the  most  complete  developments  of  the  basic 
characteristics  of  waves  in  a  free  plate  is  provided  by  Graff  based  on  work  by  Meeker 
and  Meitzler  [7].  For  simplicity  in  following  this  development,  (2.10)  through  (2.13) 
are  rewritten  in  terms  of  sines  and  cosines  (with  different  constants). 


<f>  =  (Acoskzaz  +  Bsmkzaz)e^k^  “r)  (2.14) 

=  (Ccos£zpz  +  DsmkzpZ)e*k** (2.15) 

\jry  =  (Ecoskzfiz  +  FsinJ^pZ)^**' ~“r)  (2.16) 

i|rz  =  (GcosA^z  +  HsinkzpZ)el(k^  .  (2.17) 


If  (2.2)  is  simplified  for  dependence  only  on  x  and  z  in  the  two  dimensional  problem 
under  consideration,  and  the  potentials  of  (2.14)  through  (2.17)  are  substituted  into  the 
result,  then  particle  displacements  are  given  by 

0<J>  3vji> 

“« =  a<  *  dz  ai*) 

=  [i*IOcosX’lIz.fisin*tlIz)ttlS(-£smi!8z.Fcosit|!z)]«"‘'x 
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(2.19) 


u 


y 


64;  34r 

- _  + - i 

dz  dx 

[  -  fczp  ( -  CsinAf  p  z  ^Dcosk^z) + ikx(Gcoskifiz + Hsink^z)]e  *(<vc"“,) 


= 

dz  dx  (2.20) 

=  [^a(-i4siaKzaz+jBcos^zaz)-ifcx(£cos/:zpZ+FsiiiAzpZ)]fi‘(<;xX'“,)  . 

The  basic  nature  of  wave  motion  in  elastic  solids  emerges  from  (2.18)  through  (2.20) 
in  the  decoupling  of  transverse  (SH)  motion  from  radial  and  vertical  motion  (SV/P);  uy 
depends  on  \|rx  and  yz,  while  ux  and  il,  depend  on  \\ry  and  <J). 

The  generalized  form  of  Hooke’s  law  for  the  stress  tensor,  Ty,  in  a  linearly 
elastic  solid, 


TV  Cjjpq  epq  > 


where  the  strain  tensor,  e^,  is  given  by 


ea  "  \K+ujJ>  ' 


(2.21) 


(2.22) 


reduces  in  the  homogeneous  isotropic  plate  to 


dur  duz 
x  =  (A,+2p) — -  +  X — - 
s  3z  dx 


du 


=  V- 


X  zy  =  ^ 


£  +  x 


dx 

* y 
dy 


dur) 

dz 


(2.23) 


On  the  upper  and  lower  faces  of  the  plate  (z=±h)  the  boundary  conditions 
between  plate  and  vacuum  are  simply 


19- 


(2.24) 


Equations  (2.18)  through  (2.20),  (2.23)  and  (2.24)  combine  to  form  6  equations 
in  the  eight  unknowns  A,  B,  C,  D,  E,  F,  G  and  H.  Recalling  that  Lame’s  theorem  also 
provides  that  the  divergence  of  the  vector  potential  is  zero  (2.3),  in  the  simplified 
geometry  of  the  plate  this  relation  becomes 


dx  dz 


(2.25) 


Substituting  (2.15)  and  (2.17)  into  (2.25)  yields  two  additional  equations,  for  a  total  of 
eight  equations  in  the  eight  unknowns.  In  matrix  form  this  system  of  8  equations  is 
written  as 


aC" 

0 

0 

-bsp 

fcCp 

0 

0 

A 

0 

0 

bsp 

bcp 

0 

0 

B 

0 

0 

ccp 

CSp 

0 

0 

-<fcp 

rfCp 

C 

0 

0 

ccp 

-CSp 

0 

0 

<fsp 

dct 

D 

eCa 

0 

0 

/cp 

fh 

0 

0 

E 

«C« 

0 

0 

fc  P 

~fs  p 

0 

0 

F 

0 

0 

*cp 

*5P 

0 

0 

-hsp 

*Cp 

G 

0 

0 

0 

0 

hsp 

hct 

H 

(2.26) 


where 
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a  =  (X+2p  )k2za+Xk2x 

h  =  Ki 

b  =  2i[LkJcz^ 

g  -  ikx 

C  =  4 

c„  =  cos kh 

tt  ZOL 

d  =  ikjct  p 

Sa  = 

e  =  2UcA* 

Cp  =  cosk^h 

f  = 

Sp  =  sink^h 

Using  a  number  of  straightforward  matrix  manipulations,  it  is  easy  to  show  that  the 
system  of  (2.26)  is  equivalent  to 


0Sa 

-fcSp 

0 

0 

0 

0 

0 

0 

B 

eCa 

/Cp 

0 

0 

0 

0 

0 

0 

E 

0 

0 

~dsfi 

CSp 

0 

0 

0 

0 

G 

0 

0 

gs  p 

0 

0 

0 

0 

D 

0 

0 

0 

0 

fH 

0 

0 

A 

0 

0 

0 

0 

ac« 

fcCp 

0 

0 

F 

0 

0 

0 

0 

0 

0 

CCp 

dcp 

C 

0 

0 

0 

0 

0 

0 

SCp 

/tCp 

H 

(2.28) 


A  necessary  and  sufficient  condition  for  the  existence  of  a  solution  to  the  homogeneous 
system  Ax= 0  is  that  det(A)=0  [8].  The  format  of  (2.28)  is  advantageous  because  the 
determinant  of  that  special  form  of  an  8x8  matrix  expands  to  the  simple  form 


asa  -bsp 

eca  fc  p 

where  coefficients  B  and  E  are  associated  only  with  the  elements  of  the  first,  G  and  D 
the  second,  A  and  F  the  third,  and  C  and  H  the  fourth  determinant  in  (2.29).  Thus  four 
separate  solutions  to  (2.28)  must  exist,  such  that  only  the  pair  of  coefficients  associated 


-dsp  csp 

-««  /Sp 

CCp  JCp 

-*^p  gs p 

aca  bcp 

?Cp  ACp 

=  0  , 


(2.29) 
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with  each  2x2  determinant  is  non-zero,  while  the  other  six  coefficients  vanish. 


Looking  first  at  the  coefficient  pairs  A/F  and  B/E,  if  only  A,F*0,  then 
substitution  in  (2.18)  through  (2.20)  yields 

ux  =  (ik^Acoskiaz+kzpFcoskzpz)eKk^'ut) 

uy  =  0  (2.30) 

uz  =  i-^smk^z-ik/snJc^z^^0  , 

and  if  only  B,E?^0,  the  same  substitution  produces 

ux  =  (ik^dnkzaz  p£sin^pz)e  Kk^~ut) 

uy  =  0  (2.31) 

uz  =  (KaBcoskZaz " ik^cosktpZ)e ^  . 

Solutions  (2.30)  and  (2.31)  are  polarized  in  the  longitudinal  and  vertical  directions  only; 
(2.30)  is  symmetric  with  respect  to  the  x/y  plane  while  (2.31)  is  antisymmetric  with 
respect  to  the  x/y  plane.  Additionally,  setting  the  2x2  determinant  in  (2.29)  associated 
with  each  solution’s  coefficient  pair  to  zero  generates  an  equation  which  describes  the 
required  relation  between  vertical  and  horizontal  wave  numbers  for  that  solution  to  exist, 
i.e.,  the  frequency  equation  for  that  solution.  Setting  the  2x2  determinant  for 
coefficients  A  and  F  to  zero  yields 

- esabc p -fs^aC'  =  0  ,  (2.32) 

which  can  be  rewritten  as 


gjV  =  AV-klkaK(, 

and  since 


(2.33) 
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(2.33)  becomes 


(X +2  p)£2 + kkf)  (X +2  n)(fc2  -2  \ik. 


fc2 


where 


1 


1 


~k\-2k\  kl~lkl 
p2  *  1 


,2  ,2 

^zP  kx 


(2.34) 


, 2  .2  .2 

=  *x+*zc 


and 


,2  ,2  ,2 

*P  =  kx+kzQ 


(2.35) 


tan k.nh 


4 ktkk, 


tan£  h 

ZCt 


zP"  _  ^ia  zP 

(^-^zp)2 


(236) 


(2.36)  is  th :  well  known  Rayleigh-Lamb  frequency  equation  for  symmetric  plate  waves. 
In  the  antisymmetric  case  (coefficients  B,E?*0)  the  determinant  yields  similarly 


^pft  _  (V*zp>2 


(2.37) 


Although  (2.36)  and  (2.37)  were  originally  obtained  slightly  more  than  100  years 
ago  by  Rayleigh  and  Lamb,  the  associated  frequency  spectrum  was  not  completely 
understood  until  much  later.  Figure  2-2  illustrates  the  complete  free  plate  frequency 
spectrum  as  obtained  by  Mindlin  [9]  for  a  value  of  Poisson’s  ratio  of  o=0.31.  In 
Figure  2-2  the  (positive  or  negative)  real  and  imaginary  values  of  kx=(2h/7t)kx  are 
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imaginary  _  real 

kx=(2h/7i)kx 


Figure  2-2:  Frequency  spectrum,  kx=(2h/n)kx  versus  kp=(2h/7t)kp,  for  d=0.31, 
showing  symmetric  (thick  lines)  and  anti-symmetric(dotted  lines)  modes  (from 
Mindlin  [9]). 

plotted  against  kp=(2h/rc)kp.  While  the  periodic  nature  of  the  tangent  function 
introduces  an  infinity  of  modes.  Figure  2-2  also  shows  that  except  for  the  lowest  order 
fundamental  symmetric  and  anti-symmetric  modes,  all  modes  have  a  cutoff  frequency 
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below  which  the  propagating  wavenumbers  are  imaginary  and  waves  are  evanescent  in 
the  horizontal  direction.  The  cutoff  frequencies  for  the  first  symmetric  and 
antisymmetric  modes  above  the  fundamental  ones  are  given  by 


and 


(2.38) 


respectively  [7].  Treating  a  typical  arctic  ice  sheet  (h=3m,  a=3500m/s,  P=1800m/s  [14]) 
as  a  free  plate,  below  about  300  Hz  only  the  fundamental  antisymmetric  and  symmetric 
modes  can  propagate.  Further,  for  the  case  in  which  the  wavelength  is  long  compared 
to  the  thickness  of  the  plate  (thin  plate  assumption),  i.e.,  kah— >0  and  k^h— >0,  (2.36)  can 
be  reduced  easily  to  a  simple  expression  for  the  fundamental  symmetric  mode  by 
replacing  the  tangent  functions  with  their  arguments 


(2.39) 


where  c,p  is  the  phase  speed  for  the  wave.  The  lowest  order  symmetric  mode  in  a  thin 
plate  is  called  the  longitudinal  or  plate  wave.  The  longitudinal  wave  is  non-dispersive 
in  the  thin  plate  limit  since  (2.39)  indicates  that  the  phase  speed  is  not  dependent  on 
frequency.  A  similar  expression  can  be  obtained  for  the  lowest  order  mode  of  the 
antisymmetric  case  in  the  thin  plate  limit.  By  replacing  the  tangent  function  with  its 
Taylor  series  and  discarding  higher  order  terms,  (2.37)  reduces  to 


(2.40) 


Here  cr  is  the  phase  speed  for  the  lowest  order  antisymmetric  wave,  generally  referred 


Figure  2-3:  A  symmetric  (longitudinal)  wave  ux  in  a  free  plate  seen  as  the 
superposition  of  a  pair  of  P  waves  and  a  pair  of  SV  waves  incident  on  the  faces 
of  the  plate  (after  Redwood  [11]). 

to  as  the  flexural  wave.  The  flexural  wave  is  clearly  a  dispersive  wave.  As  implied  by 
the  form  of  (2.30)  and  (2.31),  longitudinal  and  flexural  waves  in  a  free  plate  may  be 
thought  of  as  constructive  interference  of  multiply-reflected  P  and  SV  body  waves 
[10].  Figure  2-3  (after  Redwood  [11])  demonstrates  this  point  of  view  for  a 
longitudinal  wave. 

Looking  m_  it  the  remaining  two  coefficient  pairs  in  (2.28),  if  only  C,H*0,  then 
substitution  in  (2.18)  through  (2.20)  produces 

u  -  u  -  0 

1  (2  41) 

“y  =  [(k^C  +  ik^sink^e^0  , 

and  if  only  D,G/0,  then  the  same  substitution  produces 
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(2.42) 


X  Z 

uy  =  [{-k^D  +  ikxG)co^z\e^u,) 

These  two  solutions  are  transversely  polarized  SH  waves;  (2.41)  is  antisymmetric  with 
respect  to  the  x/z  plane,  and  (2.42)  symmetric  with  respect  to  the  x/z  plane.  The 
frequency  equation  from  the  2x2  determinant  in  (2.29)  associated  with  the  antisymmetric 
case  is 

ccp/tcp-dcpgcp  =  0  ,  (2.43) 

which  simplifies  to 

*, ,<&■> *,2)cosV  *°-  <2-44) 

Solutions  to  (2.44)  exist  only  for 

V  -  ■  »= 0,1  A3-  •  <2-«> 

Similarly,  for  the  symmetric  case 

kzf>(kzV  +k2x)sm2kzph  =  0  ,  (2.46) 

for  which  solutions  exist  only  for 

kzph  =  mn,  m=0,l^,3,...  .  (2.47) 

If  the  results  of  (2.45)  and  (2.47)  are  combined,  then  SH  mode  solutions  exist  for 

Kth  =  'y  •  «= 0,1,2, 3...  ,  (2.48) 

with  n  odd  the  antisymmetric  and  n  even  the  antisymmetric  modes.  (2.48)  can  be 
shifted  to  a  more  revealing  form  by  writing  it  in  terms  of  kp  and  k,. 
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kx=(2h/7c)k. 

Figure  2-4:  Frequency  spectrum,  kx  versus  k^,  for  the  symmetric(-)  and 
antisymmetric(— )  SH  modes  of  an  infinite  free  plate  (after  Graff  [7]). 


(kfh)  =  +  •  <2'49) 

(2.49)  is  plotted  in  Figure  2-4  on  axes  again  scaled  by  2h/jt.  Note  in  Figure  2-4  that  all 
but  the  fundamental  symmetric  mode  and  all  antisymmetric  modes  have  cutoff 
frequencies  given  by 


to 


P  7r  n 
2  h  1 


(2.50) 


below  which  modes  become  evanescent  (imaginary  horizontal  wavenumber);  however, 
the  fundamental  symmetric  mode,  n=0  in  (2.49),  is  a  nondispersive  propagating  wave 
independent  of  frequency.  Since  the  cutoff  frequency  given  by  (2.50)  is  relatively  high 
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for  a  typical  arctic  ice  sheet  (again,  -300  Hz  in  3  meter  ice),  it  is  likely  that  only  the 
fundamental  symmetric  SH  mode  will  be  of  concern  in  this  investigation. 


2.1.2  Elastic  Waves  in  a  Floating  Plate 

The  problem  of  a  thin  elastic  plate  over  a  liquid  half-space  differs  from  the  free 
plate  problem  only  in  the  boundary  conditions  at  z=h.  In  place  of  (2.31),  the  boundary 
conditions  are  unchanged  at  the  upper  boundary  with  vacuum. 


=  T_ 


"  V  =  °’ 


z=-h 


but  at  the  lower  boundary  they  become 


(2.51) 


T«,  =  V,  -  °- 


=  T 


zV 


Uz  =  “z,’  z=+h 
zi  h. 


(2.52) 


as  is  appropriate  at  a  liquid/solid  interface  [7].  The  subscript  1  refers  to  the  plate  and 
the  subscript  2  the  underlying  liquid  half  space  in  (2.51),  (2.52),  and  the  work  which 
follows. 

Unfortunately,  when  the  floating  plate  boundary  conditions  are  used  with 
equations  (2.18)  through  (2.20)  and  (2.23),  the  new  system  of  nine  equations  in  nine 
unknowns  cannot  be  analyzed  as  readily  as  the  free  plate  system.  The  P  and  SV  wave 
motions  in  the  floating  plate  no  longer  reduce  to  purely  symmetric  and  antisymmetric 
modes  due  to  the  presence  of  the  liquid.  Press  and  Ewing  [’2]  studied  the  P  and  SV 
modes  by  making  the  simplifying  assumption  that  X=p,  and  derived  an  exact  expression 
for  the  period  equation. 


where 


P(2Q+8coshVjJt  coshvj/O+QfisinhVjA  sinhvj/i  =  0  , 


(2.53) 
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(2.54) 


5  =  p  2g2(vi2-^)(v2-fr?) 

PlP?V2 

v  =  ,  v'  =  ikz p , 

P  =  (v/12+£^)2coshv1/i  sinhv^/z  -  4v1v/1ikf  sinhVj/i  coshv^/z , 

Q  =  (v^+fc^smliVj/i  cxyshvj/i  -  4v1Vi^coshv1/i  sinhvi/i . 

Press  and  Ewing  were  further  able  to  show  that  for  long  waves  in  a  thin  plate  (kxh 
small),  flexural  waves  analogous  to  the  free  plate’s  antisymmetric  case  exist  and 
propagate  with  a  period  equation  given  by  the  approximation 


8p  ,(M)3 

3p2 


.Jl 


2Pikxh 


(2.55) 


and  analogous  longitudinal  (symmetric  case)  waves  exist  and  propagate  with  a  period 
equation  given  by  the  approximation 


Pi  2 

2  1 — —  (1+2  i(kh)3 


(  _  ?  \  i  r  i 

MillJl  Vl?IlJl  .. 

a2P2(  cc2]  cc]{  a]) 


(2.56) 


Note  in  particular  that  the  real  part  of  die  longitudinal  wave  velocity  is  unchanged  from 
the  free  plate  (2.39),  but  that  the  liquid  half-space  adds  an  imaginary  part  which 
represents  attenuation  proportional  to  the  wavenumber  cubed.  In  the  short  wavelength 
limit  (ksh  large),  Press  and  Ewing  showed  that  (2.53)  produces  Rayleigh  waves  on  the 
free  surface  and  Rayleigh  and  Stoneley  (Scholte)  waves  at  the  ice/water  interface.  These 
interface  waves  are  not  important  for  the  ice  thicknesses  and  frequencies  of  concern  in 
this  work  and  will  not  be  discussed  further. 
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The  change  in  boundary  conditions  is  irrelevant  to  the  SH  wave,  the  derivation 
of  which  involves  neither  x27  nor  uz.  The  SH  wave  propagates  unchanged  in  a  plate 
regardless  of  the  fluid  (or  vacuum)  on  the  faces  of  the  plate.  Thus  the  solution  for  the 
SH  modes  in  the  floating  plate  is  exactly  that  for  the  free  plate  discussed  in  the  previous 
section. 


2.2  Numerical  Solutions 

Numerical  solutions  to  the  transcendental  characteristic  equations  for  the  floating 
plate  have  been  obtained  generally  by  making  simplifying  assumptions.  As  seen  in  the 
previous  section,  Press  and  Ewing  assumed  that  the  Lame  constants  were  equal  and 
looked  at  the  solution  in  the  limit  as  the  wavelength  became  very  small  or  very  large. 
A  number  of  thin  plate  theories  have  been  studied  which  model  a  fluid  loaded  plate  in 
such  a  way  as  to  account  for  only  the  lowest  flexural  and  possibly  longitudinal  modes. 
Although  Langley  [13]  has  shown  that  the  thin  plate  approximations  can  provide 
accurate  results  below  the  cutoff  frequencies  for  the  higher  order  symmetric  and 
antisymmetric  modes,  with  modem  computing  facilities  and  available  tools  these 
approximations  are  no  longer  necessary  to  achieve  quick  and  accurate  results  from  the 
exact  equations.  Stein  [14]  used  a  specialized  computational  routine  to  solve  the  P- 
SV  system  of  equations  for  the  floating  plate  numerically;  however  Schmidt  [2]  has 
developed  a  much  more  flexible  tool  to  apply  to  this  problem,  the  Seismo- Acoustic  Fast 
Field  Algorithm  for  Range-Independent  environments  (SAFARI). 

2.2.1  Full  Wavefield  Global  Matrix  Solution 

The  SAFARI  approach  to  solution  of  seismo-acoustic  propagation  in  a 
horizontally-layered  environment  is  at  its  heart  an  expansion  of  the  techniques  of  solving 
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the  depth-separated  wave  equation  originally  applied  to  acoustic  propagation  by  Pekeris 
[15]  and  later  extended  to  seismic  propagation  by  Ewing,  Jardetzky  and  Press  [10]. 
To  see  how  this  technique  works  in  SAFARI  [3],  the  wave  equation  for  compressional 
waves  (2.5)  is  rewritten  for  a  single  horizontal  layer  in  cylindrical  coordinates  with  the 
assumption  that  the  environment  is  ax i- symmetric,  range-independent  (i.e.,  two 
dimensional)  and  all  sources  are  on  the  z-axis. 


t 


V2- 


a 2  dt 2 


4>(r*l)  =  fs(z,t)  6(r) 


(2.57) 


If  the  Fourier  transform, 


F(o)  =  , 

2  nJ— 

and  the  Hankel  transform, 


(2.58) 


=  f~g(r)Jo(k/)rdr  , 


(2.59) 


are  both  applied  to  (2.57),  the  result  is  an  ordinary  differential  equation  in  z  only. 


[ d 2 


\dz: 


-(*?-*«(*))  ®(^,u)  =  - 


(2.60) 


where  lq  is  the  horizontal  wavenumber  and  ka=(£)/a  as  before.  The  solution  to  (2.60) 
is  the  depth-dependent  Green’s  function,  given  at  some  radial  frequency  0)  by 


= 

®p(^)  +  A(k^(krlz)+  \krrz). 


(2.61) 


where  Op  is  some  particular  solution  to  (2.60),  d>'  and  d>+  are  two  independent  solutions 
to  the  homogeneous  form  of  (2.60),  and  A'  and  A+  are  coefficients  determined  by  the 
boundary  conditions.  If  the  z-dependence  of  ka  is  restricted  to  cases  for  which  (2.60) 
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can  be  solved  analytically,  then  solution  of  (2.60)  in  a  multi-layer  system  amounts  to 
solving  for  the  arbitrary  coefficients  of  (2.61).  If  a  is  taken  as  constant  and  there  are 
no  sources,  for  instance,  then  solving  for  the  depth-dependent  Green’s  function  yields 

tf>(fcrTz)  =  +.A+el*“z  ,  (2.62) 

and  the  potential  is 


4>(r,z)  =  e  +  A+e*~z\jQ(krr)krdkr  .  (2.63) 

For  solid  media,  SAFARI  takes  P  as  constant  (in  each  layer)  and  solves  for 
^(r, z)  in  a  manner  similar  to  that  above;  although  as  implemented  in  SAFARI  the 
vector  \|r(r,z)  is  replaced  by  (in  cylindrical  coordinates)  \/(r,z)=9\j/0/8r.  The  equations 
for  the  displacements  and  stresses  become 


and 


«r(^)  =  -f-4 >(r»z)  + 

dr  drdz 

uz(r*)  =  ~  -irrir1*'/(r’z)  ’ 

dz  r  dr  dr 


(2.64) 


'JrA  = 


xn^ 


duT  dur 
(X+2|i)-i  + 

dz  dr 
dur  du 
p(—  +  — i)  . 

az  dr 


(2.65) 


Thus  the  depth-dependent  Green’s  function  is 


WXkrz)  =  , 


(2.66) 


and  the  appropriate  solution  for  the  potential  is 
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Vta)  =  f~  [fi  e  ~**z  +  B  teik*z]lQ(krr)dkr  .  (2.67) 

Substituting  (2.63)  and  (2.67)  into  (2.64)  and  62.65)  produces  a  set  of  four 
equations  in  the  four  unknown  coefficients  for  the  displacements, 

ufa)  =  e  -  A 

+  tt*(B  *~**Z  ~  B'e'+jf&fikjk,  (2  68) 

ufa)  =  e -  4  V*“*) 

+  +  B^^|/o(V)^r  > 

and  the  stresses. 


=  p  |o*|(2^-^)(i4  +i4  +ett“*) 

-2 ikJt^B  e  **  -Be*  J Jo(k/)kldkr  (2.69) 

=  F /J^AaO4  -  A 

+(2 k2r-k£){B-e-**z  +  B+e^)]/,(V)^r  • 

This  set  of  equations  must  be  satisfied  on  both  sides  of  each  boundary  between  layers 
such  that  the  boundary  conditions  appropriate  to  the  adjacent  layers 
(liquid/solid/vacuum)  are  satisfied.  Since  the  boundary  conditions  must  be  satisfied  at 
all  ranges,  the  kernels  of  the  integrals  in  (2.68)  and  (2.69)  must  also  satisfy  the 
boundary  conditions,  leading  immediately  to  a  linear  system  of  equations  in  the 
unknown  coefficients.  A',  A+,  B-  and  B+,  at  every  horizontal  wave  number  lq. 

In  SAFARI  the  linear  system  of  coefficients  is  solved  numerically  to  determine 
the  depth-dependent  Green’s  functions  at  desired  depths  for  a  discrete  set  of 
wavenumbers;  the  set  of  wavenumbers  chosen  must  be  sufficient  to  allow  the  numerical 
determination  of  the  inverse  Hankel  transforms  in  (2.68)  and  (2.69)  for  the  desired 
ranges  from  source  to  receiver.  The  resulting  transfer  functions  are  determined  at  a 
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discrete  set  of  frequencies,  and  again  the  set  of  frequencies  chosen  must  be  sufficient 
to  support  the  numerical  calculation  of  the  inverse  Fourier  transform  to  determine  the 
total  response  in  the  time  domain.  In  principle,  the  SAFARI  solution  is  exact,  limited 
only  by  the  accuracy  of  the  numerical  methods  used  to  solve  the  linear  systems  and 
determine  the  inverse  wavenumber  and  frequency  transforms,  and  the  need  to  describe 
the  environment  by  discrete  layers  within  which  the  wave  equation  is  separable. 
Clearly,  SAFARI  must  be  applied  with  care  and  knowledge,  in  particular  weighing  the 
requirements  of  wavenumber  and  frequency  sampling  against  the  limitations  of  machine 
memory  and  processing  time.  For  an  investigator  with  a  firm  background  in  the 
fundamentals  of  wave  propagation  and  numerical  analysis,  and  with  a  little  experience 
with  the  code,  SAFARI  is  an  ideal  tool  to  apply  to  the  problem  of  wave  propagation  in 
an  elastic  plate. 

2.2.2  Attenuation 

The  treatment  of  elastic  waves  in  a  plate  to  this  point  has  assumed  frictionless 
propagation,  but  in  reality  some  energy  is  lost  due  to  internal  friction  with  each  cycle 
of  stress.  As  developed  in  Aki  and  Richards  [4],  a  dimensionless  measure  of  this 
friction  is 


<?(«) 


a E 
2t zE 


(2.70) 


where  E  is  the  peak  strain  energy  and  -aE  is  the  energy  lost  in  each  cycle.  If 
attenuation  is  assumed  to  be  a  linear  phenomenon,  wave  amplitude.  A,  is  proportional 
to  Ve  in  a  linearly  elastic  solid.  If  Q»1  is  also  assumed,  (2.70)  can  be  rewritten 


1  _  a4 

(?(g>)  TZ  A 


(2.71) 
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i 


For  attenuation  of  a  wave  propagating  in  the  x  direction 


and  (2.71)  becomes 

dA  _  nA  _  1  kjA 

~dX  ~  <?(w)A.  ~  2  C>(o)  ’ 

for  which  the  solution  is 


(2.72) 


(2.73) 


v4(x)  =  AQe 


(2.74) 


This  relation  can  be  viewed  as  treating  the  wavenumber  of  a  propagating  wave  as  a 
complex  value,  i.e.. 


(2.75) 


The  explicit  assumptions  of  linearity  used  in  the  above  development  are  made 
in  the  SAFARI  code,  and  this  approach  to  attenuation  is  adopted  in  this  paper.  It  is 
common  in  ocean  acoustics  to  express  attenuation  in  dBA,  so  that  linear  frequency- 
dependent  attenuation  is  given  by  the  parameter  y,  where 


207iloge 

This  convention  is  used  in  SAFARI  and  throughout  this  paper. 


(2.76) 
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Chapter  3 

Experimental  Measurements 

Before  embarking  upon  detailed  discussions  of  source  location  and  inversion  in 
later  chapters,  it  is  useful  first  to  review  the  nature  of  the  experimental  measurements, 
both  to  understand  the  motivation  for  the  chosen  approach  to  propagation  analysis 
(including  the  need  for  extensive  work  with  localization  of  experimental  sources)  and 
to  observe  certain  aspects  of  the  measured  wave  forms  which  are  in  themselves 
astonishing.  This  chapter  will  first  review  the  seismo-acoustic  propagation  experiments 
which  generated  the  data  used  in  later  chapters  and  then  will  take  a  close  look  at  the 
data  collected,  including  the  time  series  of  longitudinal,  flexural  and  SH  mode  waves 
generated  by  under-ice,  underwater  explosive  charges. 

3.1  The  Experiments 

All  data  sets  included  in  this  investigation  were  collected  during  experiments 
conducted  during  March  and  April  1987  at  an  ice  camp  designated  PRUDEX,  located 
in  the  Beaufort  Sea  approximately  100  nautical  miles  north  of  Prudhoe  Bay,  Alaska. 
The  PRUDEX  ice  camp  was  established  as  a  joint  effort  of  the  Woods  Hole 
Oceanographic  Institution,  the  Massachusetts  Institute  of  Technology  and  the  Polar 
Science  Center  of  the  University  of  Washington  Applied  Physics  Laboratory  with  the 
primary  objective  of  testing  an  Arctic  Remote  Autonomous  Measurement  Platform 
(ARAMP)  data  buoy  [16];  however,  other  experiments  were  included  in  the  program. 
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Figure  3-1:  PRUDEX  hydrophone  and  geophone  array  layout  on  the  x/y  plane 
used  throughout  this  paper,  hydrophones  H0-H15  suspended  at  a  60  meter  water 
depth,  geophones  G1-G4  frozen  into  the  ice. 

Originally  intended  only  as  engineering  experiments  designed  to  verify  techniques  and 
procedures  for  full  scale  experiments  in  later  years,  seismo-acoustic  ice  propagation 
experiments  were  conducted  to  observe  the  response  of  the  PRUDEX  ice  canopy  to 
underwater  detonations.  Preliminary  review  of  the  propagation  data  sets  recorded  during 
PRUDEX  indicated  that  the  data  sets  were  of  such  high  quality  as  to  warrant  a  much 
more  detailed  investigation  than  had  been  planned  when  the  sets  were  obtained. 

Two  co-located  receiving  arrays  were  used  to  record  seismo-acoustic  data  at  the 
PRUDEX  camp.  As  shown  in  Figure  3-1,  waterborne  waves  were  received  at  a  sixteen 
element  hydrophone  array.  Each  hydrophone  was  suspended  at  a  depth  of  60  meters  at 
the  locations  shown.  No  exact  hydrophone  position  data  were  available  during  the  data 
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recordings  of  interest.  Hydrophone  response  was  sampled  at  1000  Hz  and  recorded  in 
digital  form  on  an  optical  disk  (17).  Waves  propagating  in  the  ice  were  received  by 
an  array  of  four  3-axis  geophones.  Geophones  were  frozen  in  the  ice  at  the  locations 
shown  in  Figure  3-1.  The  geophone  data  were  recorded  on  wideband  magnetic  tape  and 
later  sampled  at  1000  Hz  and  digitized  by  the  author  using  the  Woods  Hole 
Oceanographic  Institution’s  MIZEX  Analog  to  Digital  Converter  [18]. 

Seismo-acoustic  propagation  data  sets  were  generated  during  PRUDEX  with  a 
series  of  eight  sets  of  underwater  explosive  detonations  (shots),  each  set  conducted  in 
the  early  morning  (Universal  Time)  on  successive  days  in  late  March  and  early  April. 
Each  set  of  shots  was  conducted  at  a  different  location  at  ranges  of  from  300  to  nearly 
600  meters  from  the  receiving  arrays.  Each  set  consisted  generally  of  six  separate  shots 
conducted  at  source  depths  of  2,  4,  6,  8,  16  and  32  feet  below  the  ice  surface  using 
various  amounts  of  explosive  charge.  Table  3-1  summarizes  the  experimental  shots 
conducted  at  the  PRUDEX  ice  camp  for  which  data  were  available.  Of  all  shots 
completed,  two  shots  in  each  of  four  data  sets  proved  useful  for  the  propagation  analysis 
of  this  paper.  Useful  shots  were  primarily  those  which  excited  strong  flexural  and 
longitudinal  waves  in  the  ice,  as  well  as  producing  a  measurable  geophone  response  to 
the  waterborne  acoustic  wave  generated  by  the  detonation  as  the  wave  passed  the  ice 
beneath  the  geophone.  The  flexural  waves  were  essential  to  the  propagation  analysis 
and  inversion,  while  the  longitudinal  waves  and  the  waterborne  acoustic  wave  responses 
proved  to  be  necessary  for  determining  shot  location.  All  of  the  required  responses  were 
generated  by  detonations  at  depths  of  8  and  16  feet,  at  ranges  of  500-600  meters,  and 
using  1  or  2  feet  of  detonating  cord  (primacord)  as  the  explosive  source. 

Because  the  propagation  data  sets  were  planned  only  to  verify  equipment  and 
procedures  and  not  with  subsequent  analysis  in  mind,  much  of  the  supporting 
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date 

shot 

series 

shot  dcsig. 

time  (U.T.) 

depth  (feet) 

charge 

3/26/87 

B 

B2 

0537:41 

* 

6  drams 

B3 

0542:03 

* 

1  ft  cord 

B4 

0546:14 

* 

2  ft  cord 

3/27/87 

C 

Cl 

0418:37 

* 

4.5  drams 

C2 

0419:58 

* 

6  drams 

C3 

0423:29 

* 

1  ft  cord 

C4 

0427:24 

* 

2  ft  cord 

C5 

0436:00 

* 

4.5  drams 

C6 

0438:51 

* 

6  drams 

Cl 

0442:36 

♦ 

1  ft  cord 

C8 

0445:33 

* 

2  ft  cord 

C9 

0456:34 

* 

4.5  drams 

CIO 

0458:06 

♦ 

6  drams 

Cll 

0501:54 

* 

1  ft  cord 

C12 

0504:29 

* 

2  ft  cord 

3/29/87 

D 

D1 

0534:56 

* 

4.5  drams 

D2 

0537:40 

* 

6  drams 

D3 

0543:07 

* 

1  ft  cord 

D4 

0547:14 

* 

2  ft  cord 

3/31/87 

F 

F3 

0546:21 

* 

1  ft  cord 

F4 

0550:48 

* 

2  ft  cord 

4/1/87 

G 

G1 

0442:06 

2 

2  ft  cord 

G2 

0445:26 

4 

2  ft  cord 

G3 

0450:45 

8 

2  ft  cord 

G4 

0454:46 

16 

2  ft  cord 

4/2/87 

H 

HI 

0548:25 

64 

2  ft  cord 

H2 

0553:28 

32 

2  ft  cord 

H3 

0557:14 

16 

2  ft  cord 

H4 

0601:15 

8  !  2  ft  cord 

H5 

0604:39 

4 

2  ft  cord 

H6 

0606:44 

2 

2  ft  cord 

4/3/87 

I 

11 

0558:05 

64 

2  ft  cord 

12 

0602:32 

32 

2  ft  cord 

13 

0607:33 

16 

2  ft  cord 

14 

06’.  1:42 

8 

2  ft  cord 

15 

0614:32 

4 

2  ft  cord 

16 

0617:12 

2 

2  ft  cord 

Table  3-1:  Summary  of  experimental  shots  recorded  during  PRUDEX  propagation 
experiments.  (*)  indicates  no  depth  recorded  on  shot  log,  "cord"  tefers  to  primacord, 
"drams"  to  explosive  weight. 
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information  which  could  have  gready  simplified  analysis  was  not  recorded.  Principal 
among  the  deficiencies  in  recorded  data  was  that,  except  in  only  the  most  general  terms, 
neither  shot  location  nor  shot  time  was  recorded  relative  to  the  receiving  arrays.  As  a 
result,  a  mandatory  and  non-trivial  preliminary  to  detailed  propagation  analysis  was 
locating  each  shot  using  only  the  seismo-acoustic  data  received  at  the  arrays.  A  second 
deficiency  was  the  failure  to  conduct  a  geophone  calibration  either  before,  during  or 
after  the  experiment.  Additionally,  because  the  data  from  the  two  different  arrays  was 
recorded  on  completely  different  systems  and  no  attempt  was  made  to  align  the  two 
systems  precisely,  the  hydrophone  time  series  could  not  be  synchronized  with  the 
geophone  time  series  before  analysis  began. 

Further  complicating  the  propagation  analysis  were  known  discontinuities  in  the 
ice  canopy  at  PRUDEX;  the  explosive  shots  were  made  beneath  the  ice  camp’s  runway, 
a  relatively  thin  floe  of  new  ice,  while  the  receiving  arrays  were  located  on  or  below  an 
abutting  floe  of  thicker  multi-year  ice.  Again,  as  detailed  analysis  was  not  anticipated, 
the  surface  geometry  of  the  floe  abutment  was  not  surveyed,  nor  were  any  useful  ice 
thickness  or  under-ice  surveys  conducted.  In  fact,  the  only  available  ice  thickness  data 
were  recollections  of  the  personnel  who  drilled  ice  holes  in  support  of  hydrophone  array 
and  under-ice  shot  placement.  The  runway  was  informally  reported  as  being  about  one 
meter  thick,  and  the  array  ice  reported  as  being  somewhat  variable,  about  three  meters 
in  thickness. 

3.2  The  Observations 

The  purpose  of  this  section  is  to  review  the  time  series  observed  during  the 
PRUDEX  propagation  experiments  and  associate  observed  wave  forms  with  the 
appropriate  wave  types  discussed  in  Chapter  2.  Additional  characteristics  of  the 
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Figure  3-2:  Waterborne  acoustic  waves  as  recorded  at  the  output  of  PRUDEX 
array  hydrophones  (from  top  to  bottom)  #8,  #3  and  #0  in  response  to 
experimental  under-ice  shot  number  F3. 

observed  time  series  which  greatly  simplify  analysis  and  inversion  are  identified  and 
explained. 

3.2.1  Hydrophone  Data 

Figure  3-2  shows  a  typical  time  series  resulting  from  the  detonation  of  1  foot  of 
primacord  under  the  ice  at  a  water  depth  of  8  feet  and  a  range  from  array  center  of  570 
meters,  as  recorded  at  the  output  of  several  PRUDEX  array  hydrophones.  A  striking 
feature  of  each  of  the  traces  in  Figure  3-2  is  the  characteristic  shock  wave  and  bubble 
pulse  pressure  signature  of  an  underwater  explosive  detonation.  Bubble  pulses,  the 
pressure  peaks  which  follow  the  initial  detonation  shock  wave,  result  from  successive 
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oscillations  of  the  globular  mass  of  gaseous  materials  that  remains  after  the  detonation 
is  completed,  each  successive  pulse  being  weaker  than  the  preceding  one  as  remaining 
energy  is  dissipated.  Generally,  only  the  first  several  bubble  pulses  are  stiong  enough 
to  be  observable.  Although  the  peak  pressure  of  the  first  bubble  pulse  is  about  40%  that 
of  the  shock  wave,  at  lower  frequencies  the  energy  density  is  actually  higher  in  the  first 
bubble  pulse  than  in  the  shock  wave.  The  principal  peak  in  the  combined  spectrum  of 
an  underwater  detonation  occurs  at  a  frequency  of  1/T,  with  T  being  the  interval 
between  shock  wave  and  first  bubble  pulse  [19].  As  can  be  seen  in  Figure  3-2,  after 
low  pass  filtering  to  prevent  aliasing  in  the  data  acquisition  system,  the  first  bubble 
pulse  is  actually  higher  in  amplitude  than  the  shock  wave  due  to  its  higher  energy 
density  at  low  frequencies.  This  characteristic  is  common  to  most  shots  analyzed. 

3.2.2  Geophone  Data 

Geophones  in  the  PRUDEX  array  were  implanted  in  the  ice  with  principal  axes 


Figure  3-3:  PRUDEX  geophone  array  layout  showing  alignment  of  principal 
axes  on  each  geophone. 
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aligned  as  shown  in  Figure  3-3;  however,  a  much  more  useful  aspect  is  obtained  by 
resolving  the  output  of  the  x  and  y  axis  geophones  to  axes  corresponding  to  the  radial 
and  transverse  directions  relative  to  the  propagation  path  of  the  appropriate  shot. 
Figure  3-4  shows  the  time  series  corresponding  to  the  same  shot  shown  in  Figure  3-2 
as  recorded  at  one  3-axis  geophone  and  resolved  to  the  direction  of  propagation. 

The  first  arrival  on  the  radial  particle  velocity  trace  in  Figure  3-4  is  the 
longitudinal  plate  wave,  the  pulse  corresponding  to  the  detonation  shock  wave  arriving 
first  followed  by  a  larger  one  corresponding  to  the  first  bubble  pulse.  As  the 
longitudinal  wave  is  only  slightly  dispersive,  the  characteristics  of  the  underwater 
detonation  are  retained  in  its  time  series;  unfortunately  later  multiple  arrivals  generated 
due  to  complex  interactions  with  discontinuities  in  the  ice  tend  to  confuse  the  pattern 
for  the  subsequent  weaker  bubble  pulses.  The  longitudinal  wave  is  also  visible  in  the 
vertical  geophone,  although  much  less  so  than  in  the  radial  geophone,  and  later  multiple 
arrivals  dominate  the  time  series  until  the  flexural  wave  begins. 

The  flexural  wave  is  the  strongly  dispersive  wave  beginning  about  0.7  seconds 
after  detonation  in  Figure  3-4  and  continuing  to  the  end  of  the  trace.  The  flexural  wave 
dominates  the  vertical  geophone  output,  and  is  clearly  visible  in  the  radial  geophone 
output.  The  dispersive  nature  of  the  flexural  wave  destroys  the  shock  wave/bubble  pulse 
characteristic  of  the  response;  however  comparison  with  shots  made  at  the  same  location 
but  at  different  depths,  hence  with  different  bubble  pulse  intervals,  shows  that  the 
flexural  wave  behaves  as  though  its  time  of  origin  is  the  time  of  the  first  bubble  pulse, 
not  the  initial  shock.  This  correction  is  significant  when  determining  the  phase  and 
group  velocities  of  the  flexural  wave  for  comparison  with  theoretical  results. 

A  second  distinct  set  of  pulses  on  the  vertical  geophone,  arriving  at  about  0.4 
seconds  after  detonation  in  Figure  3-4,  correspond  to  the  response  of  the  ice  plate  by  the 


-44- 


time  after  detonation  (sec) 


waterborne  acoustic  pulse  arriving  at  the  underside  of  the  floe  directly  below  the 
geophones.  Since  the  travel  path  in  the  ice  is  short,  the  pulses  replicate  nearly  exactly 
the  characteristics  of  the  hydrophone  arrivals  in  Figure  3-2.  At  the  liquid/solid  interface 
the  energy  transmission  is  only  in  the  direction  of  the  normal  to  the  surface,  and  the 
waterborne  pulse  is  seen  principally  in  the  vertical  geophone;  although  discontinuities 
in  the  underside  of  the  ice  produce  a  very  small  response  in  radial  and  transverse 
geophones  as  well.  The  responses  to  the  waterborne  pulse  carry  no  information  useful 
directly  in  the  inversion  for  the  elastic  parameters  of  the  ice;  nonetheless,  it  will  be 
shown  in  Chapter  4  that  these  pulses  are  indispensable  to  accurately  determining  the 
location  of  the  experimental  shots  relative  to  the  array. 

Figure  3-5  contains  two  expanded  plots  for  the  vertical  axis  output  of  one 
geophone,  with  the  response  from  two  separate  shots  made  at  the  same  location  overlaid 
for  comparison.  One  shot  was  made  using  1  foot  of  detonating  cord  at  a  depth  of  eight 
feet,  and  the  other  shot  using  2  feet  of  cord  at  the  same  depth.  The  first  plot  in 
Figure  3-5  shows  in  detail  the  geophone  response  to  the  waterborne  wave  as  it  passes 
beneath  the  ice  on  which  the  geophone  rests,  as  well  as  the  change  in  the  bubble  pulse 
interval  due  to  the  change  in  the  size  of  the  charge.  The  second  plot  is  of  the  flexural 
wave  and  demonstrates  both  the  repeatability  of  the  wave  and  the  slight  offset  in  the 
onset  of  the  flexural  wave  introduced  by  the  difference  in  the  bubble  pulse  intervals. 

Contrary  to  all  expectations  based  on  the  theory  of  plate  wave  propagation,  the 
largest  amplitude  response  to  an  underwater  detonation  in  each  of  the  3-axis  geophones 
occurs  only  on  the  transverse  geophone  at  a  velocity  corresponding  roughly  to  the 
expected  shear  velocity  in  the  ice.  Further,  the  wave  form  of  this  transverse  response 
retains  generally  the  shock  wave/bubble  pulse  characteristic  of  the  source,  indicating  that 
this  arrival  is  at  most  only  slightly  dispersive.  The  inescapable  conclusion  to  be  drawn 
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Figure  3-5:  Vertical  geophone  #4  response  for  experimental  shots  F3(--)  and 
F4(-);  top,  the  response  to  the  waterborne  acoustic  wave  as  it  passes  under  the 
ice,  and  bottom,  the  flexural  wave. 


from  this  combination  of  characteristics  is  that  the  largest  amplitude  wave  propagating 
in  the  ice  as  a  result  of  an  under-ice,  underwater  detonation  is  the  fundamental 
symmetric  SH  mode!  This  result  is  surprising  because,  as  discussed  in  Chapter  2,  there 
should  be  no  interaction  between  SH  waves  in  the  plate  and  acoustic  waves  in  the  water, 
thus,  an  underwater  detonation  should  not  excite  SH  waves.  Additionally,  in  a 
homogeneous,  isotropic  plate  there  is  no  mechanism  to  convert  longitudinal  and  flexural 
(SV  and  P)  waves  into  SH  waves.  Either  the  detonation  does  excite  SH  waves  in  a 
manner  not  understood,  or  out-of-plane  scattering  at  discontinuities  in  the  ice  sheet 
provides  a  significant  coupling  between  SH  and  SV/P  modes.  Either  conclusion  would 
be  significant,  in  that  most  work  in  the  arctic  environment  to  date  has  ignored  the  SH 
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Figure  3-6:  Hodograph  for  the  x  and  y  axes  of  geophone  #3  showing  the 
response  to  shot  F4  at  times  from  0.1  to  0.323  seconds(-)  and  0.323  to 
0.5  seconds(-*)  after  detonation. 

mode  in  the  ice  as  insignificant  relative  to  ocean  acoustics.  The  assumption  that  the  SH 
mode  may  be  safely  discarded  appears  to  be  questionable  based  upon  these  observations. 

Figure  3-6  shows  a  typical  hodograph  constructed  from  x  and  y  geophone  outputs 
during  an  under-ice  shot.  The  hodograph  further  illustrates  the  clear  radial  polarization 
of  the  longitudinal  wave  arrival,  followed  by  the  equally  clear  transverse  polarization 
of  the  later  SH  wave  arrival.  Similar  geophone  hodographs  have  been  constructed 
during  preliminary  analysis  of  under-ice  experimental  detonations  conducted  in  1989 
during  the  MIT/WHOI  CEAREX  ice  camp  [20],  indicating  that  the  appearance  of  SH 
waves  in  the  PRUDEX  data  may  not  be  an  isolated  phenomenon. 
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Chapter  4 

Source  Location 


As  discussed  in  Chapter  3,  a  critical  part  of  the  analysis  of  propagation  data 
obtained  at  the  PRUDEX  ice  camp  is  the  two-dimensional  localization  of  the  under-ice 
explosive  shots  which  excite  seismo-acoustic  propagation  in  the  ice  cover.  Although  the 
depth  for  these  shots  is  known,  only  the  most  general  information  on  their  range  and 
bearing  from  array  center  is  available  from  records  of  the  experiment. 

Prior  to  commencing  the  analysis,  it  was  believed  that  the  information  available 
from  up  to  sixteen  hydrophones  could  provide  the  basis  for  source  localization,  and  that 
geophones  at  only  four  locations  could  contribute  little  of  additional  use.  This 
predisposition  to  rely  on  the  hydrophone  data  proved  to  be  entirely  erroneous,  not  only 
because  analysis  using  the  hydrophone  data  was  insufficiently  accurate,  but  also  because 
it  concealed  information  crucial  to  a  thorough  understanding  of  propagation  in  the  ice. 

This  chapter  reviews  the  data  available  to  the  localization  effort  and  the  routines 
used  in  localization,  then  examines  the  results  for  both  hydrophone  and  geophone  data. 
A  detailed  analysis  of  apparent  anomalies  observed  in  the  geophone  results  leads  to 
some  unexpected  conclusions  about  what  these  anomalies  reveal  about  propagation  in 
the  PRUDEX  ice  cover.  Finally,  possible  locations  for  the  source  of  the  SH  waves 
noted  in  Chapter  3  are  reviewed. 

4.1  Localization  Data 
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To  take  full  advantage  of  the  unique  shock  wave/bubble  pulse  character  of  both 
the  received  hydrophone  and  geophone  data,  the  occurrence  of  each  pulse  series  is 
identified  in  each  received  time  series  and  associated  with  the  appropriate  non-dispersive 
wave  type,  i.e.,  the  longitudinal  wave,  the  SH  wave  or  the  waterborne  pulse  response 
in  the  geophone  data,  or  the  waterborne  wave  in  the  hydrophone  data.  Times  of  each 
pulse  maximum  and  minimum  are  interpolated  to  the  nearest  0.1  msec.  Thus,  for  each 
geophone  time  series  the  basic  localization  data  consist  of  times  of  received  maxima  and 
minima  at  each  of  four  geophone  locations  for  the  shock  wave  and  the  first  several 
bubble  pulses  associated  with  each  of  three  separate  wave  types.  For  each  hydrophone 
time  series,  the  localization  data  consist  of  the  times  of  received  maxima  and  minima 
for  the  shock  wave  and  the  first  several  bubble  pulses  at  thirteen  hydrophone  locations 
(hydrophones  6,  9  and  12  in  Figure  3-1  were  not  connected  to  the  recording  system). 

The  accuracy  of  the  localization  data  so  developed  is  effected  by  the  limitations 
of  the  recording  system  and  the  interpolation  routine,  as  well  as  by  the  presence, 
principally  in  the  geophone  data,  of  interfering  waves.  In  order  to  estimate  the 
uncertainty  in  the  localization  data,  the  hydrophone  maxima  are  averaged  to  determine 
the  mean  time  difference  between  a  shot’s  shock  wave  and  each  of  its  bubble  pulses. 
Since  this  bubble  pulse  At  should  be  some  constant  value  regardless  of  wave  or  receiver 
type,  deviations  from  the  mean  value  are  used  to  develop  the  statistics  of  both 
hydrophone  and  geophone  data  so  that  an  appropriate  weight  can  be  assigned  in  the 
localization  routine. 

A  second  key  element  of  localization  data  is  the  sound  speed  in  the  first  few 
meters  of  sea  water  direcdy  beneath  the  ice.  Fortunately,  this  data  is  readily  available 
[16]  for  times  within  one  hour  of  the  shot  times  of  interest.  Figure  4-1  is  a  typical  plot 
of  temperature,  salinity,  and  sound  speed  for  the  PRUDEX  ice  camp. 
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Figure  4-1:  Temperature/Salinity/Sound  Speed  profiles  at  the  PRUDEX  ice 
camp,  31  March  1987,  0601  U.T.  (from  McPhee  [16]). 

Although  not  a  major  factor  in  the  localization,  the  depth  of  each  shot  was 
recorded  during  testing.  Informal  discussions  with  personnel  involved  in  the  experiment 
indicate  that  these  depth  measurements  were  obtained  by  lowering  the  charges  into  the 
water  on  a  marked  line  and  are  reliable  values. 

4.2  Localization  Routine 

To  provide  a  flexible  tool  with  which  to  localize  each  shot  in  the  specialized  ice 
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geometry,  a  FORTRAN  computer  routine  centered  around  a  singular  value 
decomposition  (SVD)  routine  adapted  from  Press,  et.  al.[ 21],  is  employed.  SVD 
analysis  provides  a  consistent  method  of  combining  data  points  from  many  different 
array  elements,  shots,  and  wave  forms,  while  maintaining  the  ability  to  supply 
supporting  parameters  which  might  be  known  separately,  such  as  wave  speed  and  shot 
time,  and  the  ability  to  specify  the  uncertainty  in  all  data  points  and  parameters 
individually.  A  system  of  equations  is  assembled  in  matrix  form  based  on  a  simple 
linear  relation  for  each  pulse  arrival  at  each  array  element, 

)-c  <4.i) 

where  tj‘  is  the  shot  time  for  the  jth  shot,  tjawr  the  received  time  for  a  specific  pulse  of 
wave  w  at  array  element  a  from  the  jth  shot,  ra  is  the  assumed  range  from  the  shot  to 
array  element  a  (treated  as  a  known  value),  and  l/vw  is  the  inverse  speed  of  wave  w. 
For  example,  a  simple  matrix  system  consisting  of  the  shock  waves  from  two  different 
wave  types  ( e.g .,  the  longitudinal  wave  and  the  SH  wave),  arriving  from  two  separate 
shots  conducted  at  a  single  location,  and  received  at  a  two  element  array,  becomes 
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If  one  of  the  speeds  or  shot  times  is  known,  an  additional  row  is  appended  to  matrix  A 
and  vector  b  with  this  information.  The  system  is  then  weighted  by  dividing  each  row 
of  matrix  A  and  vector  b  by  the  uncertainty  of  that  row  (making  the  reasonable 
assumption  that  the  uncertainty  in  each  equation  is  uncorrelated).  For  a  given  shot 
location,  singular  value  decomposition  solution  of  (4.2)  after  weighting  yields  the  best 
fit  to  the  data  for  shot  time  and  wave  speed.  The  uncertainty  of  each  estimate  is  also 
available  [21]. 

In  order  to  locate  a  shot  in  space  and  time,  the  above  SVD  routine  is  employed 
in  a  matched  field  approach  to  search  through  two  dimensional  space  for  a  best  fit  for 
shot  times  and  wave  velocities,  as  indicated  by  the  lowest  value  of  chi-square  (%2),  the 
sum  of  the  square  of  the  difference  between  the  modelled  array  times  and  the  data.  The 
best  fit  corresponds  to  the  most  likely  location  for  the  shot  To  evaluate  the  uncertainty 
in  the  best  fit  location,  values  of  x2  arc  determined  for  each  point  in  space  searched  by 
the  localization  routine,  converted  to  dB,  normalized  to  zero  dB  at  the  lowest  value  of 
X2  (i  e.,  converted  to  curves  of  a%2  from  the  best  fit),  and  plotted.  The  best  fit  location 
and  the  estimated  uncertainties  (assumed  to  be  Gaussian)  are  then  used  to  compute  a 
large  number  of  Monte  Carlo  simulations  of  the  data  [21].  These  Monte  Carlo 
simulations  are  supplied  to  the  localization  routine  and  their  best  fit  locations 
determined.  Simultaneously  plotting  these  best  fit  Monte  Carlo-derived  locations  with 
the  ax2  contours  resulting  from  the  original  best  fit  allows  a  straightforward  assignment 
of  confidence  limits  to  specific  dB  values  of  ax2. 

4.3  Localization  with  Hydrophone  Data 
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Although  use  of  all  available  elements  in  the  hydrophone  array  was  intended, 
when  legalization  using  hydrophone  data  was  initially  attempted,  it  became  immediately 
apparent  that  positions  recorded  for  the  elements  located  100  meters  and  further  from 
the  array  center  were  very  inaccurate.  Localization  failed  to  converge  meaningfully 
when  data  from  these  outer  hydrophones  were  incorporated,  and  attempts  to  refine  outer 
hydrophone  positions  with  data  from  different  shot  locations  on  successive  days  failed 
due  the  small  angular  separation  and  long  time  interval  between  shot  series.  Limited 
thus  to  the  inner  9  hydrophones  (of  which  8  were  connected  to  the  recording  system) 
the  localization  routine  proved  very  successful  in  refining  shot  bearing,  but  much  less 
so  in  refilling  shot  range. 

Figure  4-2  is  a  plot  of  a%2  contours  for  a  best  fit  shot  location  determined  using 
hydrophone  data,  along  with  the  best  fit  locations  for  80  Monte  Carlo  simulations  of  that 
data.  The  1  dB  ellipse,  corresponding  to  about  a  90%  confidence  limit,  has  a  major  axis 
more  than  900  meters  long.  With  an  estimated  range  from  shot  to  array  center  of  much 
less  than  this  distance,  wave  velocities  determined  using  travel  times  from  the  best  fit 
shot  location  cannot  be  specified  even  to  within  ±50%.  As  will  be  seen  in  the  next 
chapter,  inversion  of  ice  propagation  data  is  heavily  dependent  on  measured  wave 
velocities.  Clearly,  hydrophone-based  shot  locations  do  not  provide  the  accuracy  which 
reliable  inversion  requires. 

4.4  Shot  Location  using  Geophone  Data 

After  having  failed  to  determine  a  reliable  location  and  time  for  any  experimental 
shot  using  the  hydrophone  data,  localization  was  begun  using  geophone  data.  With  the 
geophone  data  it  was  hoped  that  the  information  carried  in  the  multiple  waves  could 
offset  the  limitations  of  an  array  comprised  of  only  four  more  closely  spaced  element 
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Figure  4-2:  ax2  Contours  for  the  best  fit  location  for  the  F  series  of  shots 
calculated  using  hydrophone  array  data,  plotted  with  the  best  fit  locations(+)  for 
Monte  Carlo  simulations  of  that  data. 

positions.  The  results  eventually  much  more  than  justified  the  expectations,  but  not  until 
a  thorough  analysis  and  explanation  of  some  apparent  anomalies  was  completed. 

4.4.1  Variations  in  Ice  Thickness  at  the  Receiving  Array 

The  first  anomaly  is  relatively  easy  to  understand  and  eliminate.  When  only  the 
arrivals  of  the  response  to  the  waterborne  pulse  are  processed  in  the  localization  routine, 
the  shot  bearings  correspond  well  with  those  produced  by  the  hydrophone  data,  as  is 
expected,  but  wave  velocities  are  consistently  reported  as  10-15  m/s  higher  than  the 
known  value  of  about  1435  m/s.  Since  the  relatively  low  uncertainty  of  the  pulse 
measurements  supports  more  accurate  velocity  determinations,  some  variation  in  the 
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travel  path  is  postulated.  Detailed  numerical  studies  show  that  this  problem  is  very 
likely  due  to  variations  in  the  ice  thickness  below  each  geophone  -  changes  in  pulse 
travel  time  of  as  little  as  0.3  ms,  as  would  be  caused  by  a  difference  of  1  meter  in  ice 
thickness  beneath  different  geophones,  account  for  the  velocity  difference  without 
significantly  altering  the  bearing  reported  by  the  localization  routine.  To  estimate  the 
thickness  at  each  geophone,  ice  thicknesses  at  the  geophones  are  added  as  variables  to 
the  SVD  system  in  the  localization  routine,  so  that  (4.1)  becomes 
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where  a  is  the  (known)  compressional  wave  velocity  and  2h,  the  unknown  thickness  at 
geophone  a.  To  anchor  the  estimates,  since  the  travel  times  are  only  weakly  dependent 
upon  ice  thickness,  thickness  variables  are  also  added  as  parameters  and  assigned  the 
same  average  value  and  expected  error.  With  this  modification,  the  localization  routine 
is  applied  to  the  existing  best  fit  locadon  and  values  for  the  geophone  ice  thicknesses 
determined.  These  values  are  then  used  to  correct  the  received  times,  and  the  basic 
localization  routine  is  used  to  determine  a  new  best  fit  location.  This  location  can  be 
used  recursively  to  determine  new  estimates  of  ice  thickness  and  improve  the  best 
position;  however,  the  process  generally  converges  very  rapidly.  Final  values  for  ice 
thickness  at  each  geophone  are  shown  in  Table  4-1.  Because  the  system  is  very 
sensitive  to  differences  in  thickness,  but  not  to  absolute  thickness,  only  the  differences 
from  some  average  value  are  shown. 


4.4.2  Evidence  for  Refracted  Waves 
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Geophone 

Deviation  from 

Nr 

average  thickness 

1 

+0.08  m 

2 

-0.18  m 

3 

-0.83  m 

4 

+0.97  m 

Table  4-1:  Deviations  from  the  average  ice  thickness  determined  at  each 
geophone  in  the  PRUDEX  array. 

The  basic  version  of  the  shot  location  routine  utilizes  only  the  two  most  clearly 
defined  and  essentially  non-dispersive  wave  pulse  arrivals  at  the  geophone  array,  the 
longitudinal  wave  and  the  response  to  the  impact  of  the  waterborne  pulse  on  the  ice 
under  the  geophone.  Inherent  in  the  basic  location  routine  is  the  assumption  that  the 
water  and  the  ice  plate  are  homogenous  media,  such  that  non-dispersive  waves  travel 
in  straight  lines  at  constant  speed.  As  seen  in  Figure  4-1,  this  assumption  is  valid  for 
the  waterborne  arrivals  (except  as  noted  above).  For  longitudinal  waves  traveling  in  the 
ice,  however,  the  difficulties  introduced  by  variations  in  ice  thicknesses  at  the  array 
foreshadow  the  final  conclusion  that  the  PRUDEX  ice  cover  cannot  be  treated  as  a 
single  homogeneous  plate.  When  both  waterborne  and  longitudinal  waves  are  processed 
in  the  location  routine,  the  resulting  values  are  much  higher  than  the  apparent 
measurement  uncertainty  supports;  however,  most  of  this  error  lies  in  very  poorly 
predicted  arrival  times  for  the  longitudinal  wave.  Detailed  review  of  solutions  for  all 
eight  shots  from  four  different  locations  discloses  that  in  all  cases  the  longitudinal  wave 
discrepancy  arises  because  the  longitudinal  waves  are  arriving  from  an  apparent  bearing 
about  8  to  10  degrees  to  the  right  of  the  waterborne  pulses.  A  further  review  for  any 
single  shot  location  of  the  arrivals  from  all  three  wave  types  propagating  in  the  ice 
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shows  that  each  wave  type  is  coming  from  a  different  apparent  bearing  than  the 
waterborne  pulse.  Figure  4-3  is  a  plot  of  the  apparent  direction  of  wave  arrival 
computed  for  each  wave  type  from  a  single  shot  location.  Assuming  that  all  waves 
originate  from  the  same  point.  Figure  4-3  strongly  suggests  a  family  of  waves  refracted 
in  the  horizontal  plane. 

Experimenters  returning  from  the  PRUDEX  ice  camp  described  a  relatively 
straight  ridge  in  the  ice  cover  which  marked  the  transition  from  the  thin  runway  ice 
under  which  the  experimental  shots  were  made,  to  the  thicker  ice  on  which  the  receiving 
array  was  situated.  This  description,  combined  with  the  indications  of  refraction  noted 
above,  prompted  an  investigation  ir  .o  the  possible  existence  of  the  horizontal  refraction 
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Figure  4-3:  Plot  of  PRUDEX  array  layout  showing  apparent  axis  of  arrival  of 
(1)  both  the  hydrophone  and  geophone  water  waves,  (2)  the  SH  wave,  (3)  the 
longitudinal  wave,  and  (4)  the  flexural  wave. 
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of  waves  propagating  in  the  ice  along  a  line  created  by  a  vertical  plane  separating  the 
ice  cover  into  two  homogeneous  half-plates  of  different  thicknesses. 

To  investigate  horizontal  refraction,  the  location  routine  is  further  modified  to 
determine,  for  a  given  index  of  refraction  and  ridge  line  orientation,  the  path  which  a 
wave  travels  from  source  to  receiver.  For  refracted  arrivals  (4.1)  becomes 


where  rls  is  the  distance  traveled  in  the  first  half-plate  from  the  source  to  the  ridge  line 
on  a  path  to  array  element  a,  r^  is  the  distance  traveled  in  the  second  half-plate  from 
the  ridge  line  to  array  element  a,  and  vlw  and  v2w  are  the  corresponding  wave  velocities 
in  the  two  half-plates.  Finally,  adding  an  additional  equation  to  the  system. 


(4.6) 


where  n  is  the  index  of  refraction,  establishes  the  appropriate  ratio  between  the  two 
velocities.  In  order  to  locate  th  -  most  likely  ridge  line,  a  comprehensive  search  is 
conducted  by  computing  for  each  shot  the  best  fit  shot  location  and  corresponding  value 
of  x2  for  a  broad  range  of  possible  ridge  line  orientations  and  index  of  refractions.  The 
most  likely  orientation  is  determined  by  combining  the  %2  values  for  all  shots  and 
choosing  the  orientation  that  produces  the  lowest  x2  for  all  shots  together.  Figure  4-4 
shows  the  ax2  contours  for  various  values  of  the  slope  and  y-intercept  of  the  ridge  line 
referenced  to  the  array  as  an  x/y  plane  with  the  origin  at  the  center,  and  the  positive  x 
axis  along  the  line  containing  geophone  4  and  hydrophones  4,  8  and  12.  Figure  4-5 
shows  the  orientation  of  the  best  ridge  line,  described  by  the  line 

y  -  -1.95x  +  223  meters  (4.7) 
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Figure  4-4:  a%2  contours  (in  dB)  for  the  best  fit  ridge  line  orientation, 
described  by  the  line  y=mx+b  on  an  x/y  plane  centered  on  the  horizontal  plane 
of  the  array  with  geophone  #4  on  the  x-axis. 
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Figure  4-5:  Plan  showing  on  an  x-y  plane  the  PRUDEX  geophone(*)  and 
hydrophone(+)  array,  the  best  location  for  the  F  series  of  shots(x),  the  best  fit 
ridge  line(-),  and  the  longitudinal  wave  paths  from  shot  to  each  geophone(-). 
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Figure  4-6:  Aerial  photograph  of  the  PRUDEX  iee  camp  and  array,  showing 
locations  of  identifiable  hydrophones  (geophones  and  some  hydrophones  tire  not 
visible),  the  array  axes  and  the  iee  ridge  line. 
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Figure  4-7:  a%2  contours  (in  dB)  for  the  best  fit  location  for  the  F  series  of 
shots  calculated  using  geophone  data,  plotted  with  the  best  fit  locations  for  80 
Monte  Carlo  simulations(+)  of  that  data. 

relative  to  the  array  and  the  F  series  shot  location.  For  any  given  shot  series,  the  %2 
value  for  the  best  fit  refracted  path  is  fully  6  dB  better  (lower)  than  the  %2  value  for  the 
best  fit  to  the  same  data  on  an  unrcfracted  patn.  Figure  4-6  is  the  only  available  aerial 
photograph  of  the  PRUDEX  camp  showing  both  the  ridge  line  and  the  array  layout. 
The  best  fit  ridge  line  appears  to  correspond  remarkably  well  with  the  ridge  line  in  the 
photograph. 

The  search  for  the  best  ridge  line  necessarily  includes  as  part  of  its  operation  the 
best  fit  location  for  all  shots  used  in  the  search.  Figure  4-7  is  a  plot  similar  to 
Figure  4-2  showing  the  a%2  contours  for  the  best  fit  location  for  the  F  series  shots,  as 
well  as  the  results  of  best  fit  searches  of  80  Monte  Carlo  simulations  of  the  received 


-62- 


data.  The  major  axis  of  the  90%  confidence  limit  ellipse  (corresponding  to  the  1.6  dB 
contour)  is  only  60  meters  long,  implying  an  accuracy  of  about  ±5%  in  wave  velocities 
determined  using  this  position.  In  fact,  the  slower  velocities  associated  with  the  flexural 
wave  will  tend  to  be  much  more  accurate  than  ±5%,  since  the  localization  routine 
determines  a  corresponding  shot  time  for  which  a  wave  traveling  with  the  speed  of 
sound  in  water  (1435  m/s  for  these  shots)  will  be  measured  exactly,  regardless  of  the 
error  in  range.  As  a  wave  speed  increases  or  decreases  from  this  value,  the  speed 
measurement  error  increases  accordingly.  A  wave  traveling  at  1000  m/s,  for  instance, 
will  be  measured  at  a  range  of  600  ±30  m  to  within  ±1.6%,  while  a  slower  wave 
traveling  at  500  m/s  will  be  measured  to  within  ±3%. 

Figure  4-8  provides  a  direct  comparison  of  90%  confidence  limit  ellipses  for  the 
hydrophone  and  geophone  based  best  fit  locations  for  a  shot  series.  While  a  system  to 
monitor  hydrophone  positioning,  a  more  sophisticated  processing  system,  and  a  larger 
array  can  certainly  improve  hydrophone  performance,  it  is  remarkable  that  with  all  of 
these  improvements,  the  performance  of  a  hydrophone-based  system  will  be  unlikely  to 
surpass  that  of  a  simple  system  of  four  3-axis  geophones. 

4.5  Locating  the  Source  of  the  SH  Wave 

The  SH  wave  arrivals  are  investigated  using  the  shot  location  routine  in  an 
attempt  to  determine  their  source  and  time  of  origin.  Analyzing  the  SH  wave  arrivals 
separately,  the  range  resolution  is,  as  could  be  expected  from  the  hydrophone  results, 
very  poor.  Assuming  the  point  of  origin  of  the  SH  wave  to  be  any  given  position 
between  the  shot  and  the  ridge  line,  however,  the  best  fit  time  of  SH  wave  origin 
calculated  by  the  location  routine  is  roughly  consistent  with  travel  from  the  shot  location 
to  that  point  of  origin  at  about  the  same  speed  as  the  SH  wave’s  travel  from  the  point 
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Figure  4-8:  Plots  of  the  90%  confidence  limit  ellipses  for  the  F  series  shot 
location  derived  from  geophone  and  hydrophone  data. 

of  origin  to  the  array;  i.e.,  the  SH  wave  arrivals  are  consistent  with  generation  by  the 
shot  itself,  or  with  generation  by  the  interaction  of  some  wave  traveling  at  about  the 
same  speed  as  the  SH  wave  at  some  location  near  a  line  between  the  shot  and  the  array. 

Figure  4-9  is  a  plot  of  a%2  contours  for  the  SH  wave  arrivals  determined  by 
restricting  the  time  of  origin  to  the  shot  time  determined  with  the  procedures  of  Section 
4.4.  The  intersection  of  the  1  dB  contours  for  the  two  locations  demonstrates 
qualitatively  the  compatibility  of  the  two  solutions,  and  supports  the  argument  that  the 
SH  wave  is  excited  directly  by  the  shot  itself. 

A  second  possible  source  for  the  SH  wave  is  out-of-plane  scattering  during  the 
interaction  of  the  longitudinal  wave  with  the  ridge  line.  Figure  4-10  is  a  plot  of  ax2 


-64- 


contours  for  the  SH  wave  arrivals  calculated  by  constraining  the  time  of  origin  to  the 
average  arrival  time  of  the  longitudinal  waves  at  the  ridge  line  en  route  from  the  shot 
to  the  array.  The  large  offset  of  the  center  of  the  ellipse  from  the  ridge  line  indicates 
that  the  SH  wave  arrivals  are  not  consistent  with  their  generation  during  the  interaction 
of  the  longitudinal  wave  at  the  ridge  line.  Similarly,  the  flexural  wave  interacting  at  the 
ridge  line  is  another  possible  mechanism  for  out-of-plane  scattering  and  production  of 
the  SH  wave;  however,  that  possibility  can  be  immediately  dismissed  as  the  SH  waves 
arrive  at  the  receiving  array  before  the  flexural  waves  (at  the  peak  frequency  of  20  Hz) 
arrive  at  the  ridge  line. 

A  fourth  possible  source  is  an  interaction  of  the  waterborne  wave  with  some 
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Figure  4-9:  a%2  contours  (in  dB)  generated  by  the  shot  location  routine  for  the 
SH  wave  point  of  origin(-)  assuming  the  time  of  origin  is  fixed  at  shot  time, 
with  contours  for  the  best  fit  shot  location(--)  of  section  4.4. 
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Figure  4-10:  a%2  contours  (in  dB)  generated  by  the  shot  location  ioutine  for  the 
point  of  origin  of  the  SH  wave  with  time  of  origin  fixed  at  the  average  time  of 
longitudinal  waves’  (•••)  arrival  at  the  ridge  line 

discontinuity  on  the  underside  of  the  ice.  As  the  wave  speed  in  water  is  comparable 
with  that  of  the  SH  wave,  this  possibility  cannot  be  ruled  out  immediately;  however, 
Figure  4-11,  the  a%2  contours  for  the  SH  wave  point  of  origin  assuming  that  the  SH 
wave  time  of  origin  coincides  with  the  arrival  of  the  waterborne  wave  at  the  ridge  line, 
shows  at  least  that  the  water  wave/ridge  line  interaction  is  probably  not  the  source,  again 
because  the  location  contours  for  the  time  of  that  interaction  are  offset  relatively  far 
from  the  ridge  line. 

Based  on  available  evidence,  the  source  of  the  SH  wave  cannot  be  positively 
identified.  Origin  at  or  near  the  shot  location  and  time  is  consistent  with  the  data. 
Interactions  of  the  flexural  wave,  the  longitudinal  wave  or  the  waterborne  wave  with  the 
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Figure  4-11:  a%2  contours  (in  dB)  generated  by  the  location  routine  for  the 
point  of  origin  of  the  SH  wave  with  the  time  origin  fixed  at  the  time  of  the 
waterborne  waves’(-)  arrival  at  the  ridge  line(— ). 

ridge  line  do  not  appear  to  be  likely  sources  for  the  SH  wave.  Since  the  major  known 
discontinuities  in  the  PRUDEX  ice  canopy  are  associated  with  this  ridge  line,  it  is  most 
likely  that  the  SH  wave  is  either  excited  directly  by  the  under-ice  detonation  in  an  as 
yet  unexplained  manner,  or  it  is  generated  by  out-of-plane  scattering  in  the  immediate 
vicinity  of  the  detonation,  perhaps  during  interactions  with  some  unknown  feature  on 
the  underside  of  the  ice  canopy,  in  which  case  there  is  insufficient  information  in  the 
PRUDEX  experiment  to  determine  which  wave  is  the  source. 
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Chapter  5 

Inversion  of  Propagation  Data 


Solving  the  wave  propagation  problem  in  the  manner  of  Section  2. 1 ,  computing 
particle  motion  based  on  the  equations  of  motion,  the  characteristics  of  the  material,  a 
given  geometry  and  a  known  excitation,  is  the  forward  problem  in  seismo-acoustics. 
Although  solution  of  the  forward  problem  is  seldom  easy,  the  approach  is  at  least 
straightforward  and  the  correct  solution  should  be  unique  [7],  Taking  the  measured 
particle  motions  in  response  to  a  known  source  and  processing  that  da  a  "backwards" 
through  the  appropriate  equations  to  obtain  the  unknown  elastic  and  geometric 
parameters  is  the  inverse  problem.  Inversion  of  seismo-acoustic  data  often  becomes 
more  complex  than  the  comparable  forward  problem  because  it  assumes  solution  of  the 
forward  problem  as  a  starting  point  and  must  deal  with  non-unique  solutions.  A  given 
set  of  elastic/geometric  parameters  will  produce  only  one  response  to  a  given  excitation; 
however,  it  is  very  possible  that  different  sets  of  those  parameters  will  produce  the  same 
measured  response  to  that  excitation.  The  likelihood  of  non-unique  solutions  to  the 
inversion  process  necessarily  increases  as  the  number  of  unknown  parameters  increases. 

5.1  Inversion  Parameters 

The  principal  parameters  desired  from  the  seismo-acoustic  inversion  can  be  seen 
by  inspection  of  the  equations  of  motion  in  a  linearly  elastic  solid  (2.1).  The  motion 
of  ice  particles,  hence  the  propagation  of  waves  in  the  ice,  is  dependent  upon  the  ice 
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density  p,  and  the  ice  Lame  constants  p  and  X.  Equivalently,  using  relations  (2.6)  and 
(2.7),  the  Lame  constants  may  be  expressed  as  the  compressional  and  shear  velocities 
a  and  (3. 

Although  not  expressed  directly  in  the  equations  of  motion,  internal  friction  in 
a  propagating  medium  dissipates  the  energy  of  waves  propagating  in  that  medium.  In 
most  cases  this  attenuation  must  be  known  a  priori  or  added  to  the  list  of  parameters 
to  be  determined  in  the  inversion.  In  this  work  attenuation  is  described  by  the  two 
parameters  ya  and  yp,  the  attenuation  of  compressional  and  shear  waves,  respectively, 
as  described  in  Section  2.2.2. 

Ideally,  there  are  no  unknown  geometric  factors  to  complicate  an  inversion.  In 
this  chapter  shot  location  relative  to  the  receiving  array,  as  determined  in  Chapter  4,  is 
considered  a  known  value.  The  sparse  information  available  about  ice  thickness  in  the 
vicinity  of  the  PRUDEX  ice  camp  necessitates  treating  ice  thickness  as  an  unknown  and 
including  it  in  the  inversion. 

5.2  Previous  Measurements 

Although  the  mechanical  properties  of  sea  ice  have  been  extensively  studied 
[22],  very  little  work  has  been  done  to  determine  the  low  frequency  elastic  properties 
of  the  arctic  ice  cover  in  situ.  Until  recently,  actual  measurements  have  been  limited 
to  some  early  wave  speed  measurements  in  freshwater  lake  ice  [5]  [6]  [23],  and  pack 
ice  [24|  [25|,  high  frequency  attenuation  measurements  in  glacial  ice  [26]  [27] 
and  sea  ice  [28],  wave  speed  profiling  in  both  lake  and  sea  ice  [29],  and  data 
obtained  from  small  scale  laboratory  experiments.  As  a  result,  determinations  of  the  low 
frequency  properties  of  arctic  sea  ice  were  largely  inferred  from  other  ice  environments 
or  extrapolated  from  high  frequency  laboratory  a"d  in  situ  data.  As  an  excellent 
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example  of  this  approach,  McCammon  and  McDaniel  (22]  have  employed  a 
comprehensive  summary  of  available  laboratory  and  field  measurements  to  determine 
values  for  the  attenuation  in  the  arctic  sea  ice,  for  use  in  studies  of  the  acoustic 
reflectivity  of  the  ice  cover.  Based  on  this  summary,  they  have  estimated  that 
compressional  wave  attenuation  can  be  approximated  by 

y  -  6  l0'5-a  —  (5.1) 

X 

and  (assuming  Poisson’s  ratio  to  be  constant  at  0.33)  shear  wave  attenuation  by 

Yp  =  3.610-^p  ^  .  (5.2) 

X 

Results  obtained  by  such  a  combination  of  extrapolation  and  inference  may 
accurately  represent  the  characteristics  of  seismo-acoustic  propagation  at  high 
frequencies,  but  it  is  questionable  whether  these  values  may  be  translated  to  reflect  low 
frequency  behavior  as  well.  At  low  frequencies  and  long  wavelengths,  macroscopic 
discontinuities,  such  as  cracks  and  ridges  in  an  arctic  ice  plate,  may  reduce  propagation 
speed  and  increase  attenuation  in  the  medium. 

Recently,  several  investigators  have  obtained  values  for  the  elastic  parameters  of 
arctic  ice  at  low  frequencies.  In  1986  Stein  [14][30],  as  an  adjunct  to  other  studies, 
estimated  values  for  shear  and  compressional  velocities  and  attenuations  from  earlier 
work  at  two  arctic  sites.  In  1989  Brooke  and  Ozard  [31]  completed  a  detailed  study 
of  the  elastic  properties  of  sea  ice  based  on  measurements  in  the  Slidre  Fjord  of  the 
Canadian  Archipelago  in  1986  and  1987.  The  results  of  Stein,  and  Brooke  and  Ozard 
are  summarized  in  Table  5-1. 


5.3  The  Inversion  Procedure 
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Stein 

Brooke  and  Ozard 

Smooth  Ice 

Rough  Ice 

data  date 

1980-82 

1987 

1986 

1987 

1986 

a  m/s 

3500 

NA 

NA 

NA 

NA 

c,p  m/s 

NA 

3084 

2960 

2893 

2864 

P  m/s 

1800 

1705 

1891 

1660 

1746 

ya  dB/A. 

0.46 

NA 

NA 

NA 

NA 

Yp  dB/X 

1.57 

20-40Hz 

0.32 

0.45 

2.33 

1.26 

40-80Hz 

1.00 

0.57 

2.55 

0.84 

l 

80-120Hz 

0.38 

0.49 

1.33 

0.48 

Table  5-1:  Summary  of  recent  measurements  of  the  elastic  parameters  of  arctic  sea 
ice  at  low  frequency. 


Inversion  of  the  PRUDEX  propagation  data  is  conducted  initially  with  the 
assumption  that  the  floating  ice  canopy  between  source  and  receiver  can  be  treated  as 
a  single  homogeneous  plate,  and  that  only  the  longitudinal  and  flexural  waves  are 
excited  by  the  explosive  charge.  Although  neither  of  these  assumptions  is  actually  valid, 
the  methods  of  this  procedure  serve  to  demonstrate  the  power  of  SAFARI  modeling,  and 
the  results  can  be  viewed  as  a  form  of  "average"  behavior.  The  inversion  is  then  revised 
to  reflect  the  more  complete  knowledge  of  the  plate’s  character  obtained  in  Chapter  4. 

5.3.1  Inversion  for  an  Infinite,  Homogeneous  Plate 

A  straightforward  way  to  simplify  a  seismo-acoustic  inversion  is  to  select  a  small 
subset  of  the  full  set  of  elastic  parameters  describing  the  propagating  media,  and  isolate 
for  study  a  portion  of  the  measured  response  which  is  sensitive  only  to  the  elements  of 
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this  subset.  Preliminary  study  of  the  sensitivity  of  the  flexural  wave  in  a  floating  ice 
plate  to  variations  in  the  elastic  parameters  shows  that  the  flexural  wave  is  relatively 
insensitive  to  variations  in  both  the  compressional  wave  velocity  and  the  density  of  the 
ice.  Because  of  the  insensitivity  to  the  density  of  ice,  a  nominal  value  of  .91  gm/cm3 
[32]  is  used  throughout  this  work,  and  no  attempt  is  made  to  determine  ice  density 
in  any  inversion.  More  importantly,  this  study  shows  that  the  flexural  wave  is 
dependent  only  upon  the  shear  velocity,  the  attenuation  values,  and  the  thickness  of  the 
ice  -  it  can  be  isolated  to  invert  only  for  this  more  limited  number  of  parameters.  In 
addition,  the  dispersion  curve  of  the  flexural  wave,  i.e.,  the  relation  of  the  flexural 
wave’s  group  velocity  to  its  frequency,  is  essentially  independent  of  attenuation,  and  the 
measured  dispersion  curve  can  be  inverted  for  only  the  shear  velocity  and  the  ice 
thickness. 

Adopting  an  approach  similar  to  that  used  by  Jensen  and  Schmidt  [33]  to 
determine  shear  speed  and  shear  attenuation  of  the  sea  bed  from  analogous  Scholte  wave 
data,  the  first  step  in  the  inversion  consists  of  constructing  a  dispersion  curve  for  the 
flexural  wave.  The  dispersion  curve  is  developed  using  the  flexural  wave  responses  and 
the  positions  determined  for  the  eight  experimental  shots  which  not  only  excited  a 
vigorous  flexural  wave,  but  also  excited  observable  longitudinal  waves  and  responses 
to  the  waterborne  pulse  (necessary  for  localization).  In  this  section  the  positions  used 
are  those  determined  by  treating  the  ice  sheet  as  a  homogeneous  plate.  The  dispersion 
curve  is  built  by  applying  a  moving  Fourier  transform  with  a  Hanning  window  to  the 
time  series,  with  the  necessary  window  size  determined  by  the  simple  expedient  of 
increasing  its  length  until  the  dispersion  curve  is  stabilized  [34],  All  thirty-two 
individual  curves  (eight  shots  received  at  four  vertical  geophones)  are  then  normalized 
to  a  constant  noise  value  and  combined  to  produce  a  single  best  dispersion  curve  for  the 
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assumed  homogeneous  ice  plate. 

Having  obtained  an  experimental  dispersion  curve,  the  inversion  proceeds  by 
selecting  a  likely  model  for  the  unknown  parameters,  using  the  SAFARI  numerical 
modeling  routine  to  determine  a  synthetic  dispersion  curve,  comparing  the  model  and 
the  experimental  dispersion  curves,  developing  a  correction  to  the  model,  and  recursively 
refining  the  model  until  the  model’s  curve  converges  with  the  observed  one.  If  the  ice 
thickness  is  known,  the  above  procedure  will  quickly  determine  the  correct  shear 
velocity;  unfortunately,  the  ice  thickness  at  the  PRUDEX  ice  camp  is  not  known.  For 
any  given  thickness  and  shear  velocity,  a  family  of  solutions  are  found  which  can 
reproduce  the  same  flexural  wave  velocity  at  a  given  frequency  simply  by  adjusting  the 
ice  thickness  up  or  down  and  compensating  with  an  appropriate  change  in  shear 
velocity.  While  some  difference  does  arise  between  two  such  similar  solutions  over  the 
frequency  range  of  interest  (2-60Hz),  this  difference  is  well  within  the  accuracy 
available  in  comparing  dispersion  curves.  This  problem  is  illustrated  in  Figure  5-1, 
which  shows  several  families  of  dispersion  curves  calculated  for  two  different  shear 
velocities  and  various  ice  thickness  values. 

The  uncertainty  in  shear  velocity  is  largely  eliminated  by  expanding  the  inversion 
to  include  the  longitudinal  wave.  As  the  longitudinal  wave  is  essentially  non-dispersive, 
this  expansion  necessitates  a  shift  to  the  time  domain  and  a  direct  comparison  of 
experimental  and  synthetic  time  series.  In  addition  to  providing  the  compressional  wave 
velocity,  this  expansion  has  the  added  benefit  of  allowing  the  estimation  of  the 
compressional  and  shear  wave  attenuation  as  well.  The  longitudinal  wave  is  largely 
insensitive  to  ice  thickness,  casting  some  doubt  on  its  ability  to  diminish  the  uncertainty 
in  the  shear  velocity;  however,  it  is  so  sensitive  to  shear  velocity  that  only  a  very 
limited  range  of  shear  velocities  can  combine  with  a  reasonable  compressional  wave 
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Figure  5-1:  Two  sets  of  dispersion  curves  for  flexural  waves  in  ice  at  shear 
velocities  P=1600m/s(~ )  and  P=1800m/s(-),  and  (top  to  bottom  in  each  set)  ice 
thicknesses  of  1.25,  1.15  and  1.05m. 


velocity  to  match  the  observed  longitudinal  wave.  Figure  5-2  illustrates  the  dependence 
of  the  longitudinal  wave  on  shear  and  compressional  velocities,  as  well  as  demonstrating 
its  essentially  non-dispersive  nature.  Of  equal  importance  in  reducing  the  uncertainty 
in  the  uniqueness  of  the  inversion,  matching  the  observed  flexural  wave  in  the  time 
domain  is  more  sensitive  to  errors  in  flexural  velocity  over  a  wide  frequency  range  than 
matching  the  calculated  and  observed  dispersion  curves. 

A  complication  introduced  by  the  decision  to  shift  to  comparing  synthetic  and 
observed  time  series  is  the  necessity  to  provide  an  accurate  representation  of  the 
explosive  source  used  to  excite  the  observed  waves.  A  computer  routine  based  on 
equations  provided  by  Wakeley  (35]  has  proven  to  be  very  successful  in  reproducing 
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the  acoustic  pressure  signature  measured  at  1  meter  from  a  known  underwater  explosive 
source;  although  possibly  due  to  the  effect  of  the  extreme  cold  on  the  explosives,  the 
bubble  pulse  intervals  predicted  by  Wakeley’s  routines  are  consistently  longer  than 
observed  at  PRUDEX  for  the  same  explosive  weight  and  depth.  This  discrepancy  is 
resolved  by  reducing  either  the  explosive  weight  or  the  proportionality  constant  used  in 
the  equations  slightly  from  that  provided  by  explosive  tables  [36][37]  for  a  given 
length  of  primacord  or  dram  weight  of  explosive,  such  that  the  synthetic  and  observed 
bubble  pulse  intervals  agree.  In  this  way  the  relative  spectrum  levels  of  the  real  and 
synthetic  shots  are  identical.  It  is  likely  that  the  absolute  levels  are  also  equal,  but  there 
is  insufficient  information  available  to  verify  this  assumption.  Figure  5-3  shows  the 


Figure  5-2:  Two  sets  of  group  velocity  curves  for  longitudinal  waves  in  ice 
with  compressional  velocities  of  3500m/s(-)  and  3400m/s(— )  and  shear 
velocities  (top  to  bottom  in  each  set)  of  1800,  1700  and  1600  m/s. 
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Figure  5-3:  Synthetic  time  series  for  the  pressure  signature  at  1  meter  for  an 
explosive  charge  simulating  shot  F3;  top,  sampled  at  10  KHz,  and  bottom, 
prefiltered  and  decimated  to  1000  Hz. 


Figure  5-4:  Top,  synthetic  time  series  of  Figure  5-3  filtered  to  a  2-9C  Hz  band, 
and  bottom,  spectrum  of  filtered  time  series. 
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synthetic  pressure  signature  of  the  explosive  charge  (1  foot  of  primacord)  used  in  shot 
3  of  series  F.  Note  that  after  prefiltering  and  decimation  to  a  sample  rate  of  1000  Hz 
the  first  bubble  pulse  is  larger  in  both  amplitude  and  energy  content  than  the  initial 
shock  wave,  duplicating  the  relation  seen  in  the  experimental  pulse  trains.  One  of  the 
requirements  of  the  SAFARI  pulse  calculation  routine  is  that  to  avoid  "ringing”  in  the 
output  time  series,  the  frequency  integration  routine  must  be  truncated  where  the  source 
pulse  has  a  frequency  minimum  [3].  To  m  :c  *  this  requirement,  as  well  as  to  limit  the 
computation  required  for  the  frequency  integration,  the  source  pulse  is  digitally  filtered 
[38] [39]  to  a  2-90  Hz  band,  as  shown  in  Figure  5-4.  Limiting  the  frequency 
integration  to  this  band  has  no  effect  on  the  inversion.  The  partial  spectrum  of  a  typical 


Figure  5-5:  Spectrum  of  signal  received  on  vertical  component  of  geophone  #3 
during  experimental  shot  F3,  showing  preponderance  of  energy  in  the  2-90  Hz 
band. 
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geophone  time  series,  Figure  5-5,  shows  that  most  of  the  information  carried  in  the 
signal  resides  in  that  2-90  Hz  band. 

Figure  5-6  shows  the  measured  dispersion  contours  and  Figure  5-7  the  dispersion 
contours  determined  from  a  best  fit  synthetic  time  series,  both  '  alculated  for  the  best  fit 
shot  location  determined  by  treating  the  PRUDEX  ice  cover  as  a  single  homogeneous 
plate.  Both  figures  also  include  the  exact  dispersion  curve  calculated  for  the 
homogeneous  plate’s  best  fit  inversion  parameters.  Two  figures  which  demonstrate  the 
power  of  SAFARI  pulse  modeling  are  the  plots  of  synthetic  and  observed  geophone  time 
series  for  shot  number  3  of  series  F,  Figure  5-8  for  the  horizontal  geophone  and 
Figure  5-9  for  the  vertical  geophone.  Note  in  both  figures  that  not  only  is  the  flexural 
wave  modelled  well,  but  also  the  longitudinal  wave  and  the  response  to  the  waterborne 
acoustic  pulse.  Other  arrivals  seen  after  the  longitudinal  wave,  but  before  the  flexural 
wave  are  probably  due  to  the  inhomogeneity/anisotropy  of  the  real  ice  and  are  not 
reflected  in  SAFARI  modeling.  Also  of  interest,  the  apparent  irregular  behavior  on  the 
tail  of  the  synthetic  flexural  wave  is  introduced  when  the  air  is  modeled  with  a  realistic 
sound  speed  and  density  rather  than  treated  as  a  vacuum;  however,  no  set  of  air 
parameters  modeled  the  real  response  well,  and  the  air  is  treated  as  a  vacuum  for  the 
remainder  of  this  chapter. 

5.3.2  Inversion  of  Two  Abutting  Infinite  Half-Plates 

As  discussed  in  the  previous  chapter.  The  PRUDEX  ice  plate  is  much  more 
accurately  described  as  two  half-plates  of  different  thicknesses.  Although  a  version  of 
SAFARI  able  to  handle  some  range-dependence,  including  inclusions  in  an  ice  plate  of 
a  different  thickness  than  the  rest  of  the  plate,  is  under  development  during  the  summer 
of  1990  by  Gerstoft  and  Schmidt  [40],  it  is  not  available  as  of  this  writing.  In  order 
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frequency  (Hz) 

Figure  5-6:  Observed  contours  of  spectrum  level  (in  dB  normalized  to  0  dB 
maximum)  obtained  by  combining  data  from  8  shots  at  the  PRUDEX  ice  camp, 
with  the  dispersion  curve  (-•)  calculated  for  (}=1700m/s  and  2h=1.31. 


Figure  5-7:  Synthetic  contours  of  spectrum  level  (in  dB,  normalized  to  0  dB 
maximum)  derived  from  SAFARI  time  series  calculated  for  a=3400  m/s, 
P=17(X)  m/s,  2h= 1.31m,  with  corresponding  exact  dispersion  curve  (— ). 
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Figure  5-8:  Observed  response(-)  on  radial  geophone  #3  for  PRUDEX  shot  F3; 
SAFARI  synthetic  radial  response!-*)  for  a=3400m/s,  (3=1700m/s,  2h=1.31m, 
ya=1.0dB/X,  yp=2.99dBA. 
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Figure  5-9:  Observed  response(-)  on  vertical  geophone  #3  for  PRUDEX  shot  F3; 
SAFARI  synthetic  vertical  response!- •)  for  a=3400m/s,  P=1700m/s,  2h=  1.31m, 
ya=ldBA,  Yp=2.99dBA. 


-80- 


to  provide  a  rough  inversion  simulating  the  environment  at  the  PRUDEX  ice  camp,  a 
modification  of  the  stationary  phase  method  [4]  is  used  in  conjunction  with  range 
independent  SAFARI  solutions  to  approximate  the  wave  form  received  in  the  second 
half-plate  after  detonation  in  the  first  half-plate. 

As  seen  in  Chapter  2,  the  flexural  wave  in  an  ice  plate  is  a  strongly  dispersive 
wave  of  a  single  mode.  Following  the  development  in  Aki  and  Richards  [4],  a  wave 
packet  composed  of  a  single  mode  may  be  expressed  as 

,  (5.3) 

2n 


where  |F(co)  |  is  the  spectral  density  and  <j)((o)  the  initial  phase.  In  stationary  phase 
analysis,  the  integration  path  is  along  the  real  co  axis,  and  for  large  values  of  x  and  t  the 
integrand  tot  +kxx  oscillates  rapidly,  with  each  oscillation  tending  to  cancel  the  next  in 
the  integral.  Only  at  or  near  a  saddle  point,  given  by 

-^-(-ot+k^^O  (5.4) 

do 

will  the  phase  vary  slowly  enough  to  provide  a  significant  contribution  to  the  integral. 
(5.4)  can  be  simplified  to 


x_ 

t 


do 


(5.5) 


where  u  is  the  group  velocity.  Solving  (5.5)  yields  (Os(x,t),  the  frequency  expected  to 
dominate  at  distance  x  and  time  t.  Expanding  the  phase  -cot  +  kxx  in  a  Taylor  series 
about  the  point  co=co5,  and  neglecting  higher  order  terms,  yields 


d2k . 


ot+kj  -  -uj+kx{ o)+-—j(o-oj 

4  do 


(5.6) 
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Substituting  (5.6)  into  (5.3),  and  simplifying,  results  in 


2tt 


e 


(5.7) 


where  ±  corresponds  to  d2k1/dco2  <0  or  >0,  respectively. 

In  order  to  apply  this  method  to  the  problem  of  propagation  in  two  half-plates, 
the  stationary  phase  approach  must  be  expanded  to  account  for  two  propagating  media. 
If  the  integrand  in  (5.3)  is  modified  to  reflect  propagation  for  distances  x,  and  x2  in 
media  with  horizontal  wave  numbers  and  kx2,  then  (5.4)  becomes 

-^(-ur+*jc/x1+*x2x2)  =  0  (5.8) 

da 

or 


-4-(-atl+kxIx1)  +  -f-(-u>t2+kx2K2)  -  0 

da  da 


(5.9) 


t,  +  t2  =  t 


for  which  a  solution  is 


dk 


’xl 


da 


A  *2  _  &X2 

U  da 


(5.10) 


Expanding  the  two  relations  in  (5.9)  in  a  Taylor  series  and  adding  the  results  yields 


-  w  (fj  +  f2)  +  kxlxx  +  “ 

-  ca  s(rt  + 12)  +  kxl(a  >,  +  kx2(a  )x2 


x  d2k  x  d2k 

*\  u  Kxl  x2  u  kx2 


2  da2  2  da2 


(a-as)2 . 


(5.11) 
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Substitution  in  (5.3)  gives 


fiX.l)  = 


2iz 


(  _ 2~  .  1 

*1 

\ 

d\, 

+  X2 

d2k 

u  kx2 

du>2 

da2  t 

(5.12) 


(5.12)  does  not  involve  attenuation,  so  to  provide  an  ability  to  compare  synthetic  curves 
with  experimental  ones  more  directly,  an  additional  attenuation  factor  is  used: 

y  =  e-(#A**W**»  j  (5.13) 

where  8(1  ^=7^  p -(loge  14  Ok),  x  is  the  distance  between  the  point  at  which  the  pulse 
spectral  density  is  determined  and  the  point  for  which  the  pulse  is  being  calculated,  and 
A  and  B  are  constants  determined  empirically  by  comparison  with  known  results. 

To  demonstrate  the  potential  accuracy  of  the  stationary  phase  approach  applied 
to  the  flexural  wave  in  ice,  SAFARI  is  used  to  generate  synthetic  time  series  and  phase 
and  group  velocity  curves  for  the  response  to  an  explosive  shot  for  a  given  set  of  elastic 
parameters  in  a  single  homogeneous  infinite  floating  ice  plate  at  a  ranges  of  242  m  and 
569.5  m.  The  spectral  density  of  the  short  range  shot  and  the  phase  and  group  velocity 
curves  are  supplied  to  a  computer  routine  which  uses  the  stationary  phase  approximation 
(5.7)  to  generate  the  curve  at  the  longer  range.  Figure  5-10  shows  the  short  range  time 
series,  and  Figure  5-1 1  shows  the  longer  range  time  series  as  generated  by  SAFARI  and 
as  calculated  from  the  shorter  range  pulse  using  stationary  phase  and  the  empirical 
amplitude  attenuation  of  (5.13).  Although  the  stationary  phase  result  is  not  perfect,  it 
is  nonetheless  good  enough  to  allow  at  least  an  estimation  of  best  fit  curves. 

To  apply  modified  stationary  phase  to  the  problem  of  two  abutting  half-plates, 
a  computer  routine  is  employed  which  takes  the  SAFARI-generated  phase  and  group 


time  after  detonation  (see) 

Figure  5-1 1:  SAFARI  synthetic  time  series(-)  for  parameters  of  Figure  5-10  at 
a  range  of  569.5m,  and  time  series(--)  generated  by  applying  the  method  of 
stationary  phase  to  Figure  5-10. 
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velocities  for  the  flexural  wave  in  the  two  half-plates  along  with  the  orientation  of  the 
ridge  line  separating  the  two  halves  and  determines  as  a  function  of  frequency  the  index 
of  refraction  and  the  path  from  the  shot  location  to  a  given  geophone.  The  routine  uses 
the  calculated  path,  the  supplied  spectral  density  of  the  shot  determined  from  a  time 
series  generated  for  the  travel  distance  in  the  first  plate,  and  the  phase  and  group 
velocity  curves  to  construct  using  (5.12)  the  time  series  as  received  at  the  geophone  in 
the  second  plate.  This  synthetic  time  series  for  the  flexural  wave  can  now  be  compared 
with  the  observed  time  series  to  adjust  model  parameters  and  proceed  with  the  inversion. 

Reliance  on  the  flexural  wave  for  inversion  resurrects  the  problem  of  shear 
velocity/ice  thickness  ambiguity  discussed  earlier.  To  establish  a  reasonable  anchor  on 
shear  velocity,  the  location  routine  is  used  to  investigate  the  best  fit  SH  wave  speed  in 
the  two  plates.  Assuming  that  the  SH  wave  is  excited  directly  by  the  shot  and  there  is 
no  anisotropy,  this  value  can  be  used  directly  as  the  shear  velocity  in  the  respective 
plates.  The  best  fit  consistent  with  the  longitudinal  wave’s  index  of  refraction  of  1.12 
is  p,=1590  m/s,  and  P2=1750  m/s,  for  an  index  of  refraction  of  1.10. 

The  modified  stationary  phase  procedure  introduced  above  adds  an  additional 
ambiguity  to  the  inversion  problem  at  any  given  geophone.  A  best  fit  curve  can  be 
generated  by  any  one  of  a  family  of  solutions  whose  flexural  wave  velocity  curves  and 
resulting  index  of  refraction  combine  to  produce  the  same  radial  velocity  from  source 
to  receiver.  If  the  second  plate  is  truly  a  homogeneous  half-plate  of  constant  thickness, 
then  this  ambiguity  can  be  resolved  with  the  straightforward  but  laborious  procedure  of 
locating  the  best  fit  simultaneously  at  all  four  geophone  locations;  unfortunately, 
apparent  variations  in  the  second  plate’s  thickness  have  prevented  finding  any  set  of 
parameters  for  the  second  plate  which  provide  a  good  fit  at  all,  or  even  most  geophones. 
This  problem  is  not  surprising  in  light  of  the  significant  variation  in  ice  thickness 


calculated  at  each  geophone  in  Chapter  4,  and  displayed  in  Table  4-1.  To  determine  the 
correct  solution,  the  location  routine  is  used  to  determine  the  best  index  of  refraction  for 
the  flexural  wave,  given  the  shot  location  and  the  orientation  of  the  ridge  line.  The  best 
fit  solutions  with  indexes  of  refraction  centered  about  this  value  are  chosen,  assuming 
a  constant  thickness  in  the  first  plate,  but  an  average  thickness  varying  with  geophone 
path  in  the  second. 

Figure  5-12  through  Figure  5-15  show  the  best  fit  flexural  wave  time  series  at 
each  of  the  four  vertical  geophones,  along  with  the  experimentally  observed  time  series. 
Since  no  geophone  calibration  data  is  available,  a  best  fit  geophone  calibration  factor 
of  1()'5  m/s/volt  is  applied  to  all  four  geophone  outputs  to  allow  comparison  of 
calculated  and  observed  time  series.  Solutions  at  the  four  geophones  are  summarized 
in  Table  5-2.  Figure  4-5  shows  that  the  ray  path  from  the  shot  location  to  geophone  #1 


Geophone  Nr 

Parameter 

1 

2 

3 

4 

a,  (m/s) 

3000 

Pi  (m/s) 

1590 

Ya-Yp  (dbA) 

1.0,  2.66 

2h,  (m) 

1.18 

Oj  (m/s) 

3500 

p2  (m/s) 

1750 

Ya-Yp  (dbA) 

1.0,  2.99 

2h2  (m) 

2.40 

2.37 

2.15 

2.20 

Table  5-2:  Best  compressional/shear  velocities  and  attenuations  and  plate 
thicknesses  determined  by  treating  the  PRUDEX  ice  cover  as  two  abutting  half¬ 
plates,  with  the  shot  conducted  under  plate  1,  and  the  receiving  array  on  plate  2. 
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Figure  5-12:  Observed  flexural  wave  time  series(-)  for  PRUDEX  shot  F3  at 
vertical  geophone  #4,  and  synthetic  time  series(--)  for  shot  F3  at  geophone  #4 
developed  using  the  parameters  of  Table  5-2. 
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Figure  5-13:  Observed  flexural  wave  time  series(-)  for  PRUDEX  shot  F3  at 
vertical  geophone  #3,  and  synthetic  time  series(--)  for  shot  F3  at  geophone  #3 
developed  using  the  parameters  of  Table  5-2. 
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Figure  5-14:  Observed  flexural  wave  time  series(-)  for  PRUDEX  shot  F3  at 
vertical  geophone  #1,  and  synthetic  time  seriesf--)  for  shot  F3  at  geophone  #1 
developed  using  the  parameters  of  Table  5-2. 
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Figure  5-15:  Observed  flexural  wave  time  series(-)  for  PRUDEX  shot  F3  at 
vertical  geophone  #2,  and  synthetic  time  series(-\)  for  shot  F3  at  geophone  #2 
developed  using  the  parameters  of  Table  5-2. 
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nearly  coincides  with  the  path  to  geophone  #2,  and  that  the  path  to  geophone  #3  follows 
exactly  the  path  to  geophone  #4.  Given  that  Table  5-2  shows  that  the  average  ice 
thickness  seen  by  geophones  #2  and  U3  is  less  than  that  seen  by  geophones  #1  and  #4, 
the  observation  that  the  best  fit  thickness  at  geophones  #2  and  #3  is  less  than  that  for 
geophones  #1  and  #4  is  a  confirmation  of  at  least  the  general  validity  of  theses  results. 
Even  so,  uncertainties  inherent  in  the  determinat;on  of  the  shot  location  and  the 
orientation  of  the  ridge  line,  as  well  as  in  the  application  of  the  modified  stationary 
phase  procedure  and  the  inversion  process  itself,  all  combine  to  render  the  values  of 
Table  5-2  as  no  more  than  estimates  of  elastic  parameters  at  the  PRUDEX  ice  camp. 
While  Table  5-2  should  provide  a  fair  representation  of  that  environment,  due  to  the 
complex  interaction  of  factors  involved  in  their  derivation,  it  is  not  possible  to  assign 
definite  uncertainties  to  these  parameters. 


Chapter  6 


Conclusion 

The  final  chapter  summarizes  the  significant  results  of  Chapters  3,  4  and  5,  and 
makes  some  recommendations  for  future  work. 


6.1  Summary 

Work  with  the  propagation  data  generated  at  the  PRUDEX  ice  camp  has  yielded 
a  number  of  significant  findings  which  contribute  directly  to  the  body  of  knowledge  of 
seismo-acoustic  propagation  in  the  Arctic  Ocean.  This  work  has  also  highlighted  the 
importance  of  certain  tools  in  the  furtherance  of  that  knowledge. 

6.1.1  Elastic  Parameters  of  the  Arctic  Ice 

The  values  of  bulk  compressional  and  shear  wave  speeds  obtained  for  the  thicker 
multi-year  ice  at  the  PRUDEX  ice  camp,  3500  m/s  and  1750  m/s,  respectively,  compare 
very  well  with  similar  values  obtained  by  earlier  investigators.  It  is  also  interesting  to 
note  that  the  shear  speed  measured  in  the  annual  ice,  1590  m/s,  is  considerably  lower 
than  in  the  thicker  ice,  although  not  as  low  as  some  investigators  have  predicted  [41]. 

The  work  with  the  PRUDEX  data  vigorously  supports  the  assertion  that  useful 
values  for  the  low  frequency  elastic  parameters  of  arctic  sea  ice  cannot  be  obtained  from 
laboratory  measurements  or  extrapolations  from  related  data.  Attenuation  values  of 
about  1  dB/X  for  the  compressional  wave,  and  3  dB/X  for  the  shear  wave,  as  estimated 
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in  this  work,  are  slightly  higher  than  but  consistent  with  values  reported  by  the  two 
previous  in  situ  studies  [3 1  ][  14];  however,  these  values  are  more  than  four  time  greater 
than  the  best  available  numbers  estimated  using  laboratory  and  other  data  [22]. 

6.1.2  Propagation  Mechanisms 

Study  of  the  PRUDEX  data  has  revealed  the  presence  of  strong  horizontally 
polarized  transverse  (SH)  waves  propagating  in  the  sea  ice  canopy  as  a  result  of  small 
underwater  explosive  detonations.  Since  the  theory  of  plate  wave  propagation  has  no 
mechanism  for  coupling  SH  waves  in  a  plate  with  acoustic  waves  in  an  adjacent  liquid, 
these  waves  are  entirely  unexpected  and  as  yet  unexplained.  The  PRUDEX  data  sets 
generally  support  the  contention  that  these  waves  originate  in  the  ice  canopy  at  or  very 
near  the  time  and  horizontal  location  of  the  detonation.  The  data  sets  do  not  support  SH 
wave  generation  by  out-of-plane  scattering  during  the  interaction  of  either  longitudinal, 
flexural  or  waterborne  waves  with  the  ridge  line  identified  in  the  ice  sheet. 

This  study  also  has  included  the  first  identification  of  the  horizontal  refraction 
of  a  family  of  wave  types  propagating  in  a  sheet  of  arctic  ice.  Each  of  the  wave  types 
appears  to  obey  simple  Snell’s  law  refraction  at  the  linear  abutment  between  the  two 
half-plates  which  comprise  the  ice  canopy,  refracting  at  angles  appropriate  to  the 
different  wave  speeds  in  the  two  half-plates. 

6.1.3  Analysis  Tools 

One  of  the  most  useful  lessons  highlighted  during  analysis  of  the  PRUDEX  data 
was  the  striking  superiority  of  a  simple  system  of  four  3-axis  geophones  over  a  system 
of  nine  hydrophones  in  a  larger  array.  Not  only  was  the  geophone  array  dramatically 
superior  in  localizing  the  underwater  detonations  which  excited  elastic  waves  in  the  ice, 


-91- 


it  also  allowed  the  detection  and  study  of  wave  types  and  propagation  phenomenon  not 
visible  in  the  hydrophone  data.  The  ability  of  the  geophone  array  to  isolate  several 
different  wave  types  traveling  at  different  speeds  placed  a  much  stronger  bound  on 
possible  ranges  and  bearings  from  the  array  center  to  the  shot  location  than  did  the 
hydrophone  array’s  reception  of  the  single  waterborne  wave,  despite  the  fact  that  the 
hydrophone  array  was  larger  and  the  pulse  arrivals  could  be  measured  more  accurately 
by  the  processing  system.  The  detection  of  the  SH  waves  by  the  geophone  array,  as 
well  as  the  strong  and  clear  reception  of  both  longitudinal  waves  and  flexural  waves, 
could  not  have  been  accomplished  from  hydrophone  data.  Indeed,  tomographic  studies 
of  acoustic  propagation  under  the  arctic  ice  have  shown  that  the  characteristics  of  the 
ice  cover  appear  primarily  as  second  order  effects  ( e.g .,  beam  displacement  of  the 
reflected  waterborne  pulse)  in  the  hydrophone  data  [42][43].  Clearly,  a  geophone 
array  is  a  superior  tool  for  use  in  studies  of  seismo-acoustic  propagation  in  a  localized 
section  of  sea  ice. 

Inversion  of  the  PRUDEX  data  has  also  served  to  reemphasize  the  value  of 
SAFARI  numeric  modeling  in  seismo-acoustic  propagation  problems.  In  a  homogeneous 
plate  the  results  of  Chapter  5  indicate  that  SAFARI  is  capable  of  fully  and  accurately 
reproducing  all  of  the  elements  of  the  real  seismo-acoustic  signature:  the  longitudinal 
wave,  the  flexural  wave,  and  even  the  response  to  the  waterborne  acoustic  wave  as  it 
passes  beneath  the  geophones. 

Finally,  a  potentially  useful  tool  to  extend  two-dimensional  SAFARI  to  a  range- 
dependent  environment,  the  modified  stationary  phase  approximation  of  the  flexural 
wave  (or  any  highly  dispersive  wave),  has  been  demonstrated.  This  method  serves  as 
an  effective  if  somewhat  limited  interim  fix  for  the  solution  to  propagation  in  two 
adjacent  plates  until  such  time  as  development  of  the  next  generation  of  SAFARI-like 
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algorithms  is  complete. 


6.2  Future  Work 

The  great  variation  in  the  sea  ice  elastic  parameters  determined  in  the  vicinity 
of  the  PRUDEX  ice  camp,  as  well  as  the  general  and  temporal  variability  reported 
recently  by  Brooke  and  Ozard  [31]  and  earlier  by  Hunkins  [24],  all  indicate  that  more 
work  in  this  area  is  appropriate  in  order  to  establish  a  set  of  parameters  which 
characterize  the  arctic  environment  accurately  over  a  given  area  and  season.  This  thesis 
has  shown  that  in  a  well-surveyed  homogeneous  environment  ( i.e .,  with  all  geometric 
uncertainties  eliminated),  basic  SAFARI  modeling  of  under-ice  detonations  is  readily 
capable  of  yielding  very  accurate  determinations  of  bulk  shear  and  compressional 
velocities  and  attenuations  from  geophone  data  obtained  in  situ. 

If  obtaining  geophone  measurements  in  sufficient  number  and  at  enom.'i 
locations  to  accurately  characterize  the  arctic  environment  proves  to  be  impracth  .1, 
other  approaches,  such  as  ocean  acoustic  tomography,  may  also  be  capable  of  obtaining 
average  values  of  the  elastic  parameters  over  large  areas.  High  frequency  cross-hole 
tomography  conducted  directly  in  the  ice  [29]  can  shed  additional  light  on  the  variability 
of  the  anelastic  properties  of  sea  ice,  although  extending  such  results  to  the  low 
frequencies  of  interest  in  this  work  will  remain  a  problem.  These  approaches  certainly 
warrant  further  investigation. 

Much  more  work  remains  to  be  done  to  determine  the  mechanism  which  couples 
acoustic  waves  in  the  water  with  SH  waves  in  the  plate.  An  important  element  in  this 
work  should  be  under-ice  explosive  shots  made  with  geophone  detectors  installed  not 
only  at  a  central  array,  but  spread  in  range  along  the  propagation  path.  Additionally, 
under-ice  surveys  should  be  conducted  to  identify  discontinuities  in  the  ice  canopy  not 
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visible  on  its  upper  surface.  The  information  available  in  this  expanded  experiment  may 
well  prove  vital  to  isolating  the  interaction  which  generates  the  SH  waves. 

In  order  to  further  study  seismo-acoustic  propagation  in  arctic  ice  in  a  range- 
dependent  environment,  it  will  be  necessary  to  bring  advanced  versions  of  SAFARI,  now 
in  development  [40],  to  bear  on  the  problem.  In  this  way  experimental  measurements 
will  not  necessarily  be  limited  to  strictly  homogeneous  environments,  and  data  sets  taken 
in  complex  environments,  such  as  that  of  the  PRUDEX  ice  camp,  can  be  inverted  with 
more  confidence  and  reliability  than  is  possible  with  the  limited  tools  now  available. 
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